Abundance of Deer in Victoria

Determining abundance of deer species in Victoria using camera trap and deer sign data. Supplementary code and methods to the ARI Technical Report “Abundance of Deer in Victoria”

Justin G Cally https://justincally.github.io/blog/ (Arthur Rylah Institute)https://www.ari.vic.gov.au/
07 September, 2023

About this document

This following document has been written in RMarkdown and contains embedded analyses using R and STAN programming languages. Using this RMarkdown enables us to keep all the analyses used for the report in a single repository and file. This RMarkdown has been rendered to HTML and is available

Methods used to estimate deer density

Camera detection

The density of deer at a given site can be estimated using camera-trap distance sampling (CTDS). CTDS is a modified form of distance-sampling that allows us to infer the probability a given individual will be detected within the survey area (area in front of the camera). This detection probability is a function of the distance of the individual from the camera, whereby individuals entering the camera field of view further from the camera have a lower detection probability up to a given truncation distance from the camera where detection probability is near-zero.

An underlying assumption about CTDS is that the probability a deer will be available for detection at any given point location within the camera field of view is equal. Under this assumption, for a point transect, we take into account the total area for each distance-bin area, which increases at further distance bins. However, in this study we implement a novel method that considers group size of the detected species in the availability calculations. For larger groups, CTDS should account for the availability of the closest individual rather than the availability of all individuals. When group size increases it is more likely the the closest individual from the group is closer to the camera. If we assume that the triggering of the camera is dependent upon the closest individual than we must adjust our estimated availability to account for variable group sizes. If we do not adjust for group size and only use the distance to the closest individual for our DS models then we will likely underestimate the detection rate as distances will be smaller than reality. Alternatively, if we record distances to multiple individuals in the same photo and take an average or model them independently we would likely overestimate detection probability because individuals at further distances are recorded because a closer individual has triggered the camera trap. For more information on how group size is likely to effect CTDS abundance estimates a simulation study has been written here: https://justincally.github.io/blog/posts/2022-11-10-ctds-for-groups/

Setup

Load in relevant packages for analysis, additionally, connect to the database. Camera trap and site data is stored on the database.

Show code
library(cmdstanr)
library(dplyr)
library(sf)
library(bayesplot)
library(loo)
library(gt)
library(terra)
library(caret)
library(tidyterra)
library(tidyr)
library(VicmapR)
library(kableExtra)
library(brms)
library(Distance)
library(ggrepel)
library(pROC)
options(brms.backend = "cmdstanr")
check_cmdstan_toolchain(fix = TRUE, quiet = TRUE)
register_knitr_engine(override = FALSE)
options(mc.cores=8)

#### Database connection ####
con <- weda::weda_connect(password = keyring::key_get(service = "ari-dev-weda-psql-01",
                                                username = "psql_user"), username = "psql_user")

Custom Functions

Additional functions used in the data preparation, modelling and analysis are available in the /functions directory.

Show code
# Source internal functions
sapply(list.files("functions", full.names = T), source, verbose = F)

Prepare Data

Wrangle and format data for the STAN models for the various species.

Scope of models

Outline which species should be modeled, and which projects to source data from.

Show code
# Species to run model for.
deer_species_all <- c("Cervus unicolor", "Dama dama", "Cervus elaphus", "Axis porcinus")
species_names <- c("Sambar Deer", "Fallow Deer", "Red Deer", "Hog Deer")
# Projects to select.
project_short_name <- c("hog_deer_2023", "StatewideDeer")
# Buffer for data extraction.
spatial_buffer <- 1000
# Covariate Rasters
raster_files <- "/Volumes/DeerVic\ Photos/Processed_Rasters"
# raster_files <- "data/prediction_raster"
prediction_raster <- "data/prediction_raster/statewide_raster.tif"
# For the integrated model we place limits on the maximum density of deer to integrate over. In cases where there are no detections on the camera this is limited to 15 per km2. In cases where deer were detected on the camera this is expanded to 50. We believe this is sufficiently high. Very high values of these will be less efficient. 
n_max_no_det <- 15
n_max_det <- 30

Camera locations

Download the camera locations from the database, this table outlines the locations and the deployment history of the cameras.

Show code
cams_curated <- tbl(con, dbplyr::in_schema("camtrap", "curated_camtrap_operation")) %>%
  dplyr::filter(ProjectShortName %in% !!project_short_name) %>% # Only retrieve cam records for hog deer 2023
  dplyr::collect() %>%
  dplyr::arrange(SiteID) %>%
  sf::st_as_sf(., coords = c("Longitude", "Latitude"), crs = 4283) 

n_site <- nrow(cams_curated)

Formulas for detection and abundance

Here we outline formulas to be used in the models. The formulas account for the various fixed-effect parameters.

Show code
#### Model formulas ####

#### Transect Formula: Survey Only ####
transect_formula <- ~Survey

#### Abundance Formula Options #### 
ab_formula_1 <- ~ scale(sqrt(BareSoil)) + scale(sqrt(Nitrogen)) + scale(log(PastureDistance)) + scale(BIO15) + scale(sqrt(ForestEdge)) + 
  scale(sqrt(MRVBF)) + scale(sqrt(SLOPE)) + scale(sqrt(NonPhotosyntheticVeg))

ab_formula_2 <- ~ scale(sqrt(BareSoil)) + scale(sqrt(Nitrogen)) + scale(log(PastureDistance)) + scale(BIO15) + scale(sqrt(ForestEdge)) + scale(sqrt(MRVBF)) 

ab_formula_3 <- ~ scale(sqrt(BareSoil)) + scale(sqrt(Nitrogen)) + scale(log(PastureDistance)) + scale(BIO15) + scale(sqrt(ForestEdge)) 

ab_formula_4 <- ~ 1

#### Detection Formula: Distance-sampling ####
det_formula <- ~ scale(HerbaceousUnderstoryCover) # average detection across all sites 

Create model data

Using the prepare_model_data() function we generate the data for the STAN model. This function will:
1. Download data from the database 2. Format that data to match the distance-sampling bins
3. Organised the counts into sites and group sizes
4. Generate model matrix of the various sub-models (distance-sampling, abundance and transect) 5. Generate data for the prediction process 6. Generate data for the random effect (bioregion)
7. Generate data for regional predictions (indexing based on DEECA regions)

Show code
if(!file.exists("data/multispecies_data.rds")) {
  
crown_land <- readRDS("data/crown_land.rds")
  
multispecies_data <- prepare_model_data_multispecies(species = deer_species_all,
                                 projects = project_short_name,
                     buffer = spatial_buffer,
                     detection_formula = det_formula,
                     abundance_formula = ab_formula_1,
                     transect_formula = transect_formula,
                     con = con,
                     raster_dir = raster_files,
                     prediction_raster = prediction_raster,
                     n_max_no_det = n_max_no_det,
                     n_max_det = n_max_det,
                     crown_land = crown_land,
                     evaltransects = TRUE, 
                     filter_behaviour = TRUE)

multispecies_data$keyfun <- 1 # 0 = HN, 1 = HZ
multispecies_data$raw_data[is.na(multispecies_data$raw_data)] <- 0
multispecies_data$transects[is.na(multispecies_data$transects)] <- 0

evc_groups <- readRDS("data/evc_groupnames.rds")

multispecies_data <- c(multispecies_data, evc_groups)

saveRDS(multispecies_data, "data/multispecies_data.rds")
} else {
  multispecies_data <- readRDS("data/multispecies_data.rds")
}

Model Execution

MCMC settings

Below we list the MCMC setting for our model. We run models on eight parallel chains for 400 iterations each (200 warm-up and 200 sampling). These setting provide us with 1,600 posterior draws (8 x 200).

Show code
# STAN settings
ni <- 200 # sampling iterations
nw <- 200 # warm-up iterations 
nc <- 8 # number of chains

Read in Models

These are the models used in the analysis. The first is an integrated model that requires transect and camera data. The second is a count-only model that just requires the camera data.

Show code
functions {

    /* Half-normal function
   * Args:
   *   sigma: sigma
   *   midpoints: midpoints
   * Returns:
   *   detection probability
   */
    vector halfnorm(real sigma, vector midpoints) {
      int bins = rows(midpoints);
      vector[bins] p_raw; // detection probability

      p_raw = exp(-(midpoints)^2/(2*sigma^2));
      return p_raw;
    }

      /* Hazard function
   * Args:
   *   sigma: sigma
   *   theta: theta
   *   midpoints: midpoints
   * Returns:
   *   detection probability
   */
    vector hazard(real sigma, real theta, vector midpoints) {
      int bins = rows(midpoints);
      vector[bins] p_raw; // detection probability

      p_raw = 1 - exp(-(midpoints/sigma)^(-theta));
      return p_raw;
    }

  vector prob_dist(real sigma, real theta, int keyfun, vector midpoints){
  int bins = rows(midpoints);
  vector[bins] out; // detection probability

  if(keyfun == 0){
    out = halfnorm(sigma, midpoints);
  } else if(keyfun == 1){
    out = hazard(sigma, theta, midpoints);
  }
  return out;
  }
}
data {
  int<lower=0> N;                      // number of observations
  int<lower=0> S;                      // number of species
  real delta;                          // bin width
  int<lower=1> n_site;                 // site
  int<lower=1> n_distance_bins;        // distance bins
  int<lower=1> n_gs;                   // number of group sizes
  vector[n_gs] gs;                     // group sizes
  vector[n_distance_bins] midpts;      // midpt of each interval
  real<lower=1> max_distance;         // truncation distance
  int<lower=1> max_int_dist;          // max distance as an integer
  real<lower=0> theta_frac;           // fraction of camera view
  array[n_site] int effort;           // effort
  array[n_site, n_gs, S] int n_obs;      //number of observations
  array[n_site, n_distance_bins, n_gs] int y; // observations matrix

  // summary of whether species is known to be present at each site
  int<lower = 0, upper = 1> any_seen[S, n_site];

  // number of surveys at each site
  int<lower = 0> n_survey[n_site];

  // availability prior information
  real<lower=0> bshape;               // availability shape
  real<lower=0> bscale;               // availability scale

  // camera trap detection parameters
  int<lower=0> det_ncb;                 // number of covariates for detection model
  matrix[n_site, det_ncb] det_model_matrix; // detection model matrix
  array[n_gs, n_distance_bins] real pa ; // averaged availability for multiple individuals

  // Abundance/occupancy model matrix
  int<lower = 1> m_psi;
  matrix[n_site, m_psi] X_psi;
  // negbinom scale
  // real reciprocal_phi_scale;

  //transect level information
  int<lower=1> trans;                  // total number of transects across all sites for all methods
  array[S, trans] int<lower = 0, upper = 1> y2; // transect based binomial detections
  int<lower = 0, upper = trans> start_idx[n_site];
  int<lower = 0, upper = trans> end_idx[n_site];
  int<lower=1> trans_det_ncb;           // number of covariates for transect detection model
  matrix[trans, trans_det_ncb] trans_pred_matrix; // transect detection model matrix
  int<lower=1>  n_max[n_site, S]; // max for poisson RN

  // Prediction data
  int<lower=1> npc;                 // number of prediction grid cells
  matrix[npc, m_psi] X_pred_psi; // pred matrix
  vector[npc] prop_pred; //offset
  // bioregion RE
  int<lower=1> np_bioreg;
  int<lower=1> site_bioreg[n_site];
  int<lower=1> pred_bioreg[npc];
    // region data
  int<lower=1> np_reg;
  int<lower=1> site_reg[n_site];
  int<lower=1> pred_reg[npc];

  // key function, 0 = halfnorm
  int keyfun;
}

transformed data {
  // vector[n_distance_bins] pi; // availability for point
  vector[n_site] survey_area;
  array[S, n_site] real cam_seen;
    for(i in 1:n_site) {
      survey_area[i] = theta_frac * effort[i] * pi() * (max_distance/1000)^2;
      for(s in 1:S) {
      cam_seen[s,i] = sum(n_obs[i,,s]);
      }
    }

}

parameters {
 // abundance parameters
  //array[n_site] simplex[n_gs] eps_ngs; // random variation in group size
  array[S] vector[m_psi] beta_psi;
  vector[det_ncb] beta_det;
  real log_theta;
  // transect detection parameters
  vector[trans_det_ncb] beta_trans_det;
  // temporal availability parameters
  real<lower=0, upper=1> activ;
  // bioregion RE
  array[S] real<lower=0> bioregion_sd;
  array[S] vector[np_bioreg] bioregion_raw;
  // eps group size params
  array[S] vector[n_gs] zeta;
  array[S] matrix[n_site, n_gs] eps_raw;
  array[S] real<lower=0> grp_sd;
  // od
  array[S] real od_mu;
  real odRNmu;
}

transformed parameters {
  // random effects
  array[S] vector[np_bioreg] eps_bioregion; // bioregion random effect
  // distance parameters
  array[n_site] real log_sigma;
  array[n_site] real sigma;
  array[n_site] vector[n_distance_bins] p_raw; // detection probability
  array[n_site, n_distance_bins, n_gs] real p_raw_scale; // detection probability for point independence model and accounting for availability
  array[n_site, n_distance_bins, n_gs] real<upper=0> log_p_raw; // log detection probability accounting for availability
  array[n_site, n_gs] real log_p; // log probability with summed exponential for the multinomial logit model
  array[n_site, n_gs] real<lower=0,upper=1> p; // log probability with summed exponential for the multinomial logit model
  // abundance parameters
  array[S, n_site, n_gs] real<lower=0> lambda;
  // activity parameters
  real log_activ = log(activ);
  vector[trans] logit_trans_p = trans_pred_matrix * beta_trans_det; // observation process model
  // lp_site for RN model
  array[S, n_site] real lp_site;
  vector<lower=0,upper=1>[trans] r = inv_logit(logit_trans_p);
  array[S] vector[n_site] log_lambda_psi;
  // eps group size
  array[S, n_site] simplex[n_gs] eps_ngs; // random variation in group size
  array[S, n_site] vector[n_gs] epsi;
  array[S] matrix[n_site, n_gs] eps_site;
  real<lower=0> theta = exp(log_theta);
  array[S] real od; // bioregion random effect
  real odRN = exp(odRNmu);  // bioregion random effect
  array[S, n_site] real<lower=0, upper=1> pbar;

  for(s in 1:S) {
    for(b in 1:np_bioreg) {
    eps_bioregion[s,b] = bioregion_sd[s] * bioregion_raw[s,b];
  }
  od[s] = exp(odRNmu + od_mu[s]);
  }



for(n in 1:n_site) {
  log_sigma[n] = det_model_matrix[n,] * beta_det;
  sigma[n] = exp(log_sigma[n]);
  p_raw[n] = prob_dist(sigma[n], theta, keyfun, midpts);
  for (i in 1:n_distance_bins) {
      // assuming a half-normal detection fn from line
      // p_raw[n,i] = exp(-(midpts[i])^2/(2*sigma[n]^2)); // half-normal function (pg 419 of AHM)
      // hazard function
      // p_raw[n,i] = 1 - exp(-(midpts[i]/theta)^(-sigma[n])); // hazard function (pg 507 of AHM)
    for(j in 1:n_gs) {
      p_raw_scale[n,i,j] = p_raw[n,i]*pa[j,i]; //  pr(animal occurs and is detected in bin i)
      log_p_raw[n,i,j] = log(p_raw_scale[n,i,j]);
      }
  }

  for(s in 1:S) {
  vector[n_max[n,s]+1] Nlik;
  vector[n_max[n,s]+1] gN;
  vector[n_max[n,s]+1] pcond;
  log_lambda_psi[s,n] = X_psi[n,] * beta_psi[s] + eps_bioregion[s,site_bioreg[n]];

  for(j in 1:n_gs) {
    eps_site[s, n,j] = grp_sd[s] * eps_raw[s,n,j];
    epsi[s,n,j] = exp(zeta[s,j] + eps_site[s,n,j]);
  }
  eps_ngs[s,n,] = epsi[s,n,]/sum(epsi[s,n,]);

  // p-bar

  for(k in 1:(n_max[n,s]+1))
    Nlik[k] = exp(neg_binomial_2_log_lpmf(k-1 | log_lambda_psi[s,n], odRN));

  gN = Nlik/sum(Nlik);

  for(k in 1:(n_max[n,s]+1))
    pcond[k] = (1 - (1 - r[start_idx[n]])^(k-1)) * gN[k];

  pbar[s,n] = sum(pcond);

  }

  for(j in 1:n_gs) {
    log_p[n,j] = log_sum_exp(log_p_raw[n,,j]);
    p[n,j] = exp(log_p[n,j]);
    // model site abundance
    for(s in 1:S) {
    lambda[s,n,j] = exp(log_lambda_psi[s,n] + log_p[n,j] + log_activ + log(eps_ngs[s,n,j])) .* survey_area[n];
    }
  }

// Royle-Nichols implementation in STAN (looping over possible discrete values of N)
// https://discourse.mc-stan.org/t/royle-and-nichols/14150
// https://discourse.mc-stan.org/t/identifiability-across-levels-in-occupancy-model/5340/2
    for(s in 1:S) {
if (n_survey[n] > 0) {
  vector[n_max[n,s]] lp;
    if(any_seen[s,n] == 0){ // not seen
      lp[1] = log_sum_exp(neg_binomial_2_log_lpmf(0 | log_lambda_psi[s,n], odRN), neg_binomial_2_log_lpmf(1 | log_lambda_psi[s,n], odRN) + bernoulli_lpmf(y2[s,start_idx[n]:end_idx[n]] | r[start_idx[n]:end_idx[n]]));
    } else {
      lp[1] = neg_binomial_2_log_lpmf(1 | log_lambda_psi[s,n], odRN) + bernoulli_lpmf(y2[s,start_idx[n]:end_idx[n]] | r[start_idx[n]:end_idx[n]]);
    }
    for (k in 2:n_max[n,s]){
      lp[k] = neg_binomial_2_log_lpmf(k | log_lambda_psi[s,n], odRN) + bernoulli_lpmf(y2[s,start_idx[n]:end_idx[n]] | 1-(1-r[start_idx[n]:end_idx[n]])^k);
    }
    lp_site[s,n] = log_sum_exp(lp);
    } else {
    lp_site[s, n] = 0;
    }
  }
    }
}

model {

  for(s in 1:S) {
    beta_psi[s] ~ normal(0, 3); // prior for intercept in poisson model
    bioregion_sd[s] ~ normal(0,2);
    bioregion_raw[s,] ~ normal(0,1);
    to_vector(eps_raw[s,,]) ~ std_normal();
    grp_sd[s] ~ normal(0, 1);
    zeta[s,] ~ normal(0, 2);
    // od
    od_mu[s] ~ normal(0,1);
  }
  odRNmu ~ normal(0,1);
  beta_trans_det ~ normal(0, 2); // beta for transect detection
  beta_det ~ normal(0, 4); // prior for sigma
  activ ~ beta(bshape, bscale);  //informative prior
  log_theta ~ normal(0,2); // prior for theta

  for(n in 1:n_site) {
  for(j in 1:n_gs) {
  for(s in 1:S) {
  target += neg_binomial_2_lpmf(n_obs[n,j,s] | lambda[s,n,j], od[s]);
        }
  y[n,,j] ~ multinomial_logit(to_vector(log_p_raw[n,,j]));
  }
  for(s in 1:S) {
  target += lp_site[s, n];
  }
}
}

generated quantities {
  array[S, n_site,n_gs] real n_obs_pred;
  array[S, n_site, n_gs] real n_obs_true;
  array[S, n_site] real N_site;
  array[S, n_site] real N_site_pred;
  array[n_site, max_int_dist+1] real DetCurve;
  array[n_site, n_gs] real log_lik1;
  array[S, n_site, n_gs] real log_lik2;
  array[n_site, n_gs] real log_lik2_site;
  array[S, n_site] real log_lik2_species;
  array[n_site] real log_lik;
  array[n_site] real log_lik_det;
  //array[S, n_site] real Site_lambda;
  real av_gs[S];
  array[S] simplex[n_gs] eps_gs_ave;
  array[S, npc] real pred;
  //array[S, npc] real pred_trunc;
  array[S, np_reg] real Nhat_reg;
  array[S, np_bioreg] real Nhat_bioreg;
  // array[np_reg] real Nhat_reg_design;
  real Nhat[S];
  //real Nhat_trunc[S];
  real Nhat_sum;
  int trunc_counter[S];
  trunc_counter[S] = 0;

for(s in 1:S) {
  eps_gs_ave[s] = exp(zeta[s])/sum(exp(zeta[s]));
  av_gs[s] = sum(gs .* eps_gs_ave[s]); //average group size
}


for(n in 1:n_site) {
  for(j in 1:n_gs) {
  log_lik1[n,j] = multinomial_logit_lpmf(y[n,,j] | to_vector(log_p_raw[n,,j])); //for loo
  for(s in 1:S) {
  log_lik2[s,n,j] = neg_binomial_2_lpmf(n_obs[n,j,s] | lambda[s,n,j], od[s]); //for loo
  n_obs_true[s, n, j] = gs[j] * (neg_binomial_2_log_rng(log_lambda_psi[s,n] + log(eps_ngs[s,n,j]), od[s]));
  n_obs_pred[s, n,j] = gs[j] * neg_binomial_2_rng(lambda[s,n,j], od[s]);
  }
  log_lik2_site[n, j] = log_sum_exp(log_lik2[,n,j]);
    }
    // get loglik on a site level
    log_lik_det[n] = log_sum_exp(log_lik1[n,]);
    log_lik[n] = log_sum_exp(log_sum_exp(log_lik_det[n],
    log_sum_exp(log_lik2_site[n,])), log_sum_exp(lp_site[,n]));
      for(s in 1:S) {
    //Site_lambda[s,n] = exp(log_lambda_psi[s,n]);
    log_lik2_species[s, n] = log_sum_exp(log_sum_exp(log_lik2[s,n,]), lp_site[s,n]);
    N_site[s,n] = sum(n_obs_true[s,n,]);
    N_site_pred[s,n] = sum(n_obs_pred[s,n,]);
      }

  // loop over distance bins
  if(keyfun == 0) {
  for(j in 0:max_int_dist) { // get DS predictions for distance 0 to max bin distance
    DetCurve[n, j+1] = exp(-(j+0.5)^2/(2*sigma[n]^2)); // half normal
    }
  } else if(keyfun == 1) {
    for(j in 0:max_int_dist) { // get DS predictions for distance 0 to max bin distance
    DetCurve[n, j+1] =  1 - exp(-((j+0.5)/sigma[n])^(-theta)); //hazard rate
    }
  }
}

for(s in 1:S) {
for(i in 1:np_reg) Nhat_reg[s,i] = 0;
for(i in 1:np_bioreg) Nhat_bioreg[s,i] = 0;

for(i in 1:npc) {
  pred[s,i] = neg_binomial_2_log_rng(X_pred_psi[i,] * beta_psi[s] + eps_bioregion[s, pred_bioreg[i]], od[s]) * prop_pred[i] * av_gs[s]; //offset
  if(pred[s,i] > max(N_site[s,])) {
    trunc_counter[s] += 1;
  }
  // for Hog Deer limit to Gippsland region
  if(pred_reg[i] != 2 && s == 4) {
    pred[s,i] = 0;
  }
  Nhat_reg[s,pred_reg[i]] += pred[s,i];
  Nhat_bioreg[s,pred_bioreg[i]] += pred[s,i];
}
Nhat[s] = sum(pred[s,]);
//Nhat_trunc[s] = sum(pred_trunc[s,]);
}
Nhat_sum = sum(Nhat);
}
Show code
functions {

    /* Half-normal function
   * Args:
   *   sigma: sigma
   *   midpoints: midpoints
   * Returns:
   *   detection probability
   */
    vector halfnorm(real sigma, vector midpoints) {
      int bins = rows(midpoints);
      vector[bins] p_raw; // detection probability

      p_raw = exp(-(midpoints)^2/(2*sigma^2));
      return p_raw;
    }

      /* Hazard function
   * Args:
   *   sigma: sigma
   *   theta: theta
   *   midpoints: midpoints
   * Returns:
   *   detection probability
   */
    vector hazard(real sigma, real theta, vector midpoints) {
      int bins = rows(midpoints);
      vector[bins] p_raw; // detection probability

      p_raw = 1 - exp(-(midpoints/sigma)^(-theta));
      return p_raw;
    }

  vector prob_dist(real sigma, real theta, int keyfun, vector midpoints){
  int bins = rows(midpoints);
  vector[bins] out; // detection probability

  if(keyfun == 0){
    out = halfnorm(sigma, midpoints);
  } else if(keyfun == 1){
    out = hazard(sigma, theta, midpoints);
  }
  return out;
  }
}
data {
  int<lower=0> N;                      // number of observations
  int<lower=0> S;                      // number of species
  real delta;                          // bin width
  int<lower=1> n_site;                 // site
  int<lower=1> n_distance_bins;        // distance bins
  int<lower=1> n_gs;                   // number of group sizes
  vector[n_gs] gs;                     // group sizes
  vector[n_distance_bins] midpts;      // midpt of each interval
  real<lower=1> max_distance;         // truncation distance
  int<lower=1> max_int_dist;          // max distance as an integer
  real<lower=0> theta_frac;           // fraction of camera view
  array[n_site] int effort;           // effort
  array[n_site, n_gs, S] int n_obs;      //number of observations
  array[n_site, n_distance_bins, n_gs] int y; // observations matrix

  // summary of whether species is known to be present at each site
  int<lower = 0, upper = 1> any_seen[S, n_site];

  // number of surveys at each site
  int<lower = 0> n_survey[n_site];

  // availability prior information
  real<lower=0> bshape;               // availability shape
  real<lower=0> bscale;               // availability scale

  // camera trap detection parameters
  int<lower=0> det_ncb;                 // number of covariates for detection model
  matrix[n_site, det_ncb] det_model_matrix; // detection model matrix
  array[n_gs, n_distance_bins] real pa ; // averaged availability for multiple individuals

  // Abundance/occupancy model matrix
  int<lower = 1> m_psi;
  matrix[n_site, m_psi] X_psi;
  // negbinom scale
  // real reciprocal_phi_scale;

  //transect level information
  int<lower=1> trans;                  // total number of transects across all sites for all methods
  array[S, trans] int<lower = 0, upper = 1> y2; // transect based binomial detections
  int<lower = 0, upper = trans> start_idx[n_site];
  int<lower = 0, upper = trans> end_idx[n_site];
  int<lower=1> trans_det_ncb;           // number of covariates for transect detection model
  matrix[trans, trans_det_ncb] trans_pred_matrix; // transect detection model matrix
  int<lower=1>  n_max[n_site, S]; // max for poisson RN

  // Prediction data
  int<lower=1> npc;                 // number of prediction grid cells
  matrix[npc, m_psi] X_pred_psi; // pred matrix
  vector[npc] prop_pred; //offset
  // bioregion RE
  int<lower=1> np_bioreg;
  int<lower=1> site_bioreg[n_site];
  int<lower=1> pred_bioreg[npc];
    // region data
  int<lower=1> np_reg;
  int<lower=1> site_reg[n_site];
  int<lower=1> pred_reg[npc];
  // evc data
  int<lower=1> np_evc;
  int<lower=1> site_evc[n_site];
  int<lower=1> pred_evc[npc];


  // key function, 0 = halfnorm
  int keyfun;
}

transformed data {
  // vector[n_distance_bins] pi; // availability for point
  vector[n_site] survey_area;
  array[S, n_site] real cam_seen;
    for(i in 1:n_site) {
      survey_area[i] = theta_frac * effort[i] * pi() * (max_distance/1000)^2;
      for(s in 1:S) {
      cam_seen[s,i] = sum(n_obs[i,,s]);
      }
    }

}

parameters {
 // abundance parameters
  //array[n_site] simplex[n_gs] eps_ngs; // random variation in group size
  array[S] vector[m_psi] beta_psi;
  vector[det_ncb] beta_det;
  real log_theta;
  // transect detection parameters
  vector[trans_det_ncb] beta_trans_det;
  // temporal availability parameters
  real<lower=0, upper=1> activ;
  // bioregion RE
  array[S] real<lower=0> bioregion_sd;
  array[S] vector[np_bioreg] bioregion_raw;
  // evc RE
  real<lower=0> evc_sd;
  array[S] vector[np_evc] evc_raw;
  // eps group size params
  array[S] vector[n_gs] zeta;
  array[S] matrix[n_site, n_gs] eps_raw;
  array[S] real<lower=0> grp_sd;
  // od
  array[S] real od_mu;
  real odRNmu;
}

transformed parameters {
  // random effects
  array[S] vector[np_bioreg] eps_bioregion; // bioregion random effect
  array[S] vector[np_evc] eps_evc; // bioregion random effect
  // distance parameters
  array[n_site] real log_sigma;
  array[n_site] real sigma;
  array[n_site] vector[n_distance_bins] p_raw; // detection probability
  array[n_site, n_distance_bins, n_gs] real p_raw_scale; // detection probability for point independence model and accounting for availability
  array[n_site, n_distance_bins, n_gs] real<upper=0> log_p_raw; // log detection probability accounting for availability
  array[n_site, n_gs] real log_p; // log probability with summed exponential for the multinomial logit model
  array[n_site, n_gs] real<lower=0,upper=1> p; // log probability with summed exponential for the multinomial logit model
  // abundance parameters
  array[S, n_site, n_gs] real<lower=0> lambda;
  // activity parameters
  real log_activ = log(activ);
  vector[trans] logit_trans_p = trans_pred_matrix * beta_trans_det; // observation process model
  // lp_site for RN model
  array[S, n_site] real lp_site;
  vector<lower=0,upper=1>[trans] r = inv_logit(logit_trans_p);
  array[S] vector[n_site] log_lambda_psi;
  // eps group size
  array[S, n_site] simplex[n_gs] eps_ngs; // random variation in group size
  array[S, n_site] vector[n_gs] epsi;
  array[S] matrix[n_site, n_gs] eps_site;
  real<lower=0> theta = exp(log_theta);
  array[S] real od; // bioregion random effect
  real odRN = exp(odRNmu);  // bioregion random effect
  array[S, n_site] real<lower=0, upper=1> pbar;

  for(s in 1:S) {
    for(b in 1:np_bioreg) {
    eps_bioregion[s,b] = bioregion_sd[s] * bioregion_raw[s,b];
  }
    od[s] = exp(odRNmu + od_mu[s]);
  for(k in 1:np_evc) {
  eps_evc[s,k] = evc_sd * evc_raw[s,k];
}
  }

for(n in 1:n_site) {
  log_sigma[n] = det_model_matrix[n,] * beta_det;
  sigma[n] = exp(log_sigma[n]);
  p_raw[n] = prob_dist(sigma[n], theta, keyfun, midpts);
  for (i in 1:n_distance_bins) {
      // assuming a half-normal detection fn from line
      // p_raw[n,i] = exp(-(midpts[i])^2/(2*sigma[n]^2)); // half-normal function (pg 419 of AHM)
      // hazard function
      // p_raw[n,i] = 1 - exp(-(midpts[i]/theta)^(-sigma[n])); // hazard function (pg 507 of AHM)
    for(j in 1:n_gs) {
      p_raw_scale[n,i,j] = p_raw[n,i]*pa[j,i]; //  pr(animal occurs and is detected in bin i)
      log_p_raw[n,i,j] = log(p_raw_scale[n,i,j]);
      }
  }

  for(s in 1:S) {
  vector[n_max[n,s]+1] Nlik;
  vector[n_max[n,s]+1] gN;
  vector[n_max[n,s]+1] pcond;
  log_lambda_psi[s,n] = X_psi[n,] * beta_psi[s] + eps_bioregion[s,site_bioreg[n]] + eps_evc[s,site_evc[n]];

  for(j in 1:n_gs) {
    eps_site[s, n,j] = grp_sd[s] * eps_raw[s,n,j];
    epsi[s,n,j] = exp(zeta[s,j] + eps_site[s,n,j]);
  }

  eps_ngs[s,n,] = epsi[s,n,]/sum(epsi[s,n,]);

  // p-bar
  for(k in 1:(n_max[n,s]+1))
    Nlik[k] = exp(neg_binomial_2_log_lpmf(k-1 | log_lambda_psi[s,n], odRN));

  gN = Nlik/sum(Nlik);

  for(k in 1:(n_max[n,s]+1))
    pcond[k] = (1 - (1 - r[start_idx[n]])^(k-1)) * gN[k];

  pbar[s,n] = sum(pcond);

  }

  for(j in 1:n_gs) {
    log_p[n,j] = log_sum_exp(log_p_raw[n,,j]);
    p[n,j] = exp(log_p[n,j]);
    // model site abundance
    for(s in 1:S) {
    lambda[s,n,j] = exp(log_lambda_psi[s,n] + log_p[n,j] + log_activ + log(eps_ngs[s,n,j])) .* survey_area[n];
    }
  }

// Royle-Nichols implementation in STAN (looping over possible discrete values of N)
// https://discourse.mc-stan.org/t/royle-and-nichols/14150
// https://discourse.mc-stan.org/t/identifiability-across-levels-in-occupancy-model/5340/2
    for(s in 1:S) {
if (n_survey[n] > 0) {
  vector[n_max[n,s]] lp;
    if(any_seen[s,n] == 0){ // not seen
      lp[1] = log_sum_exp(neg_binomial_2_log_lpmf(0 | log_lambda_psi[s,n], odRN), neg_binomial_2_log_lpmf(1 | log_lambda_psi[s,n], odRN) + bernoulli_lpmf(y2[s,start_idx[n]:end_idx[n]] | r[start_idx[n]:end_idx[n]]));
    } else {
      lp[1] = neg_binomial_2_log_lpmf(1 | log_lambda_psi[s,n], odRN) + bernoulli_lpmf(y2[s,start_idx[n]:end_idx[n]] | r[start_idx[n]:end_idx[n]]);
    }
    for (k in 2:n_max[n,s]){
      lp[k] = neg_binomial_2_log_lpmf(k | log_lambda_psi[s,n], odRN) + bernoulli_lpmf(y2[s,start_idx[n]:end_idx[n]] | 1-(1-r[start_idx[n]:end_idx[n]])^k);
    }
    lp_site[s,n] = log_sum_exp(lp);
    } else {
    lp_site[s, n] = 0;
    }
  }
    }
}

model {

  for(s in 1:S) {
    beta_psi[s] ~ normal(0, 3); // prior for intercept in poisson model
    bioregion_sd[s] ~ normal(0,2);
    bioregion_raw[s,] ~ normal(0,1);
    to_vector(eps_raw[s,,]) ~ std_normal();
    grp_sd[s] ~ normal(0, 1);
    zeta[s,] ~ normal(0, 2);
    // od
    od_mu[s] ~ normal(0,1);
    // evc re
    evc_raw[s,] ~ normal(0,1);
  }
  evc_sd ~ normal(0,1);
  odRNmu ~ normal(0,1);
  beta_trans_det ~ normal(0, 2); // beta for transect detection
  beta_det ~ normal(0, 4); // prior for sigma
  activ ~ beta(bshape, bscale);  //informative prior
  log_theta ~ normal(0,2); // prior for theta

  for(n in 1:n_site) {
  for(j in 1:n_gs) {
  for(s in 1:S) {
  target += neg_binomial_2_lpmf(n_obs[n,j,s] | lambda[s,n,j], od[s]);
        }
  y[n,,j] ~ multinomial_logit(to_vector(log_p_raw[n,,j]));
  }
  for(s in 1:S) {
  target += lp_site[s, n];
  }
}
}

generated quantities {
  array[S, n_site,n_gs] real n_obs_pred;
  array[S, n_site, n_gs] real n_obs_true;
  array[S, n_site] real N_site;
  array[S, n_site] real N_site_pred;
  array[n_site, max_int_dist+1] real DetCurve;
  array[n_site, n_gs] real log_lik1;
  array[S, n_site, n_gs] real log_lik2;
  array[n_site, n_gs] real log_lik2_site;
  array[S, n_site] real log_lik2_species;
  array[n_site] real log_lik;
  array[n_site] real log_lik_det;
  //array[S, n_site] real Site_lambda;
  real av_gs[S];
  array[S] simplex[n_gs] eps_gs_ave;
  array[S, npc] real pred;
  //array[S, npc] real pred_trunc;
  array[S, np_reg] real Nhat_reg;
  // array[np_reg] real Nhat_reg_design;
  real Nhat[S];
  //real Nhat_trunc[S];
  real Nhat_sum;
  int trunc_counter[S];
  trunc_counter[S] = 0;

for(s in 1:S) {
  eps_gs_ave[s] = exp(zeta[s])/sum(exp(zeta[s]));
  av_gs[s] = sum(gs .* eps_gs_ave[s]); //average group size
}


for(n in 1:n_site) {
  for(j in 1:n_gs) {
  log_lik1[n,j] = multinomial_logit_lpmf(y[n,,j] | to_vector(log_p_raw[n,,j])); //for loo
  for(s in 1:S) {
  log_lik2[s,n,j] = neg_binomial_2_lpmf(n_obs[n,j,s] | lambda[s,n,j], od[s]); //for loo
  n_obs_true[s, n, j] = gs[j] * (neg_binomial_2_log_rng(log_lambda_psi[s,n] + log(eps_ngs[s,n,j]), od[s]));
  n_obs_pred[s, n,j] = gs[j] * neg_binomial_2_rng(lambda[s,n,j], od[s]);
  }
  log_lik2_site[n, j] = log_sum_exp(log_lik2[,n,j]);
    }
    // get loglik on a site level
    log_lik_det[n] = log_sum_exp(log_lik1[n,]);
    log_lik[n] = log_sum_exp(log_sum_exp(log_lik_det[n],
    log_sum_exp(log_lik2_site[n,])), log_sum_exp(lp_site[,n]));
      for(s in 1:S) {
    //Site_lambda[s,n] = exp(log_lambda_psi[s,n]);
    log_lik2_species[s, n] = log_sum_exp(log_sum_exp(log_lik2[s,n,]), lp_site[s,n]);
    N_site[s,n] = sum(n_obs_true[s,n,]);
    N_site_pred[s,n] = sum(n_obs_pred[s,n,]);
      }

  // loop over distance bins
  if(keyfun == 0) {
  for(j in 0:max_int_dist) { // get DS predictions for distance 0 to max bin distance
    DetCurve[n, j+1] = exp(-(j+0.5)^2/(2*sigma[n]^2)); // half normal
    }
  } else if(keyfun == 1) {
    for(j in 0:max_int_dist) { // get DS predictions for distance 0 to max bin distance
    DetCurve[n, j+1] =  1 - exp(-((j+0.5)/sigma[n])^(-theta)); //hazard rate
    }
  }
}

for(s in 1:S) {
for(i in 1:np_reg) Nhat_reg[s,i] = 0;

for(i in 1:npc) {
  pred[s,i] = neg_binomial_2_log_rng(X_pred_psi[i,] * beta_psi[s] + eps_bioregion[s, pred_bioreg[i]] + eps_evc[s, pred_evc[i]], od[s]) * prop_pred[i] * av_gs[s]; //offset
  if(pred[s,i] > max(N_site[s,])) {
    //pred_trunc[s,i] = max(N_site[s,]);
    trunc_counter[s] += 1;
  }
  // for Hog Deer limit to Gippsland region
  if(pred_reg[i] != 2 && s == 4) {
    pred[s,i] = 0;
  }
  Nhat_reg[s,pred_reg[i]] += pred[s,i];
}
Nhat[s] = sum(pred[s,]);
//Nhat_trunc[s] = sum(pred_trunc[s,]);
}
Nhat_sum = sum(Nhat);
}
Show code
model_negbin_ms <- cmdstan_model(here::here("stan", "count_det_nondet_negbin_ms.stan"))
model_negbin_ms_evc <- cmdstan_model(here::here("stan", "count_det_nondet_negbin_ms_evc.stan"))

Fit models

Detection Models

Using the distance package we compare several detection functions (hazard and half-normal) with the inclusion of a parameter (herbaceous understorey cover), which may impact detection rates. The top model (as chosen by AIC) from this selection will be included in our full STAN model.

Show code
filter_behaviour <- TRUE
snapshot_interval <- 2
  # Camera information
  theta <- 40 * pi / 180 # camera angle in radians

  ### Get camera information for the camera sites
  ### For cameras that have problems you need to use the date from when the problem started
  cams_curated2 <- dplyr::tbl(con, dbplyr::in_schema("camtrap", "curated_camtrap_operation")) %>%
    dplyr::filter(ProjectShortName %in% !!project_short_name) %>% # Only retrieve cam records for hog deer 2023
    dplyr::collect() %>%
    dplyr::mutate(DateTimeDeploy = as.POSIXct(DateTimeDeploy),
                  DateTimeRetrieve = as.POSIXct(DateTimeRetrieve),
                  Tk = as.numeric(DateTimeRetrieve - DateTimeDeploy, units = "secs"), #seconds
                  Tk_prob = dplyr::coalesce(as.numeric(as.POSIXct(Problem1_to,
                                                                  format = "%Y-%m-%d %H:%M:%OS") -
                                                         as.POSIXct(Problem1_from, format = "%Y-%m-%d %H:%M:%OS"),
                                                       units = "secs"), 0),
                  Tk_adj = Tk-Tk_prob,
                  Tkt = Tk_adj / snapshot_interval, # snapshot moments: every second second
                  Effort = (Tk_adj*theta)/(snapshot_interval * pi * 2)) %>%
    dplyr::arrange(SiteID)


  ### Get the deer records (hog deer)
  camtrap_records_deer <- dplyr::tbl(con, dbplyr::in_schema(schema = "camtrap", table = "curated_camtrap_records")) %>%
    dplyr::filter(scientific_name %in% deer_species_all & ProjectShortName %in% !!project_short_name) %>%
    dplyr::select(Species = scientific_name, SiteID, Distance = metadata_Distance, size = metadata_Multiples, Date, Time, Behaviour = metadata_Behaviour) %>%
    dplyr::collect() %>%
    dplyr::mutate(Distance = dplyr::case_when(Distance == "NA" | is.na(Distance) ~ "999",
                                              TRUE ~ Distance)) %>%
    dplyr::mutate(Time = as.POSIXct(Time, format = "%H:%M:%OS"),
                  Time_n = as.numeric(Time, units = "secs"),
                  Behaviour = na_if(x = Behaviour, y = "NA")) %>%
    dplyr::rowwise() %>%
    dplyr::mutate(DistanceMod = list(stringr::str_split(Distance, pattern = "_&_")[[1]])) %>%
    dplyr::mutate(Distance = DistanceMod[which.min(as.numeric(stringr::str_extract(DistanceMod, pattern = "[0-9]+")))]) %>%
    dplyr::mutate(Distance = dplyr::case_when(Distance == "0-2.5" ~ "0 - 2.5",
                                              Distance == "999" ~ NA_character_,
                                              TRUE ~ Distance)) %>%
    dplyr::ungroup() %>%
    dplyr::filter(Time_n %% snapshot_interval == 0) %>% #& #snapshot moment interval of 2s
    {if(filter_behaviour) dplyr::filter(., is.na(.data[["Behaviour"]])) else .} %>% # filter out behaviors such as camera or marker interaction
    dplyr::group_by(SiteID, Time_n, Species) %>%
    dplyr::slice(1) %>% # if two photos occur in one second take only one (snapshot moment = 2)
    dplyr::ungroup()

  ### Tidy format for the records
  dcount<- camtrap_records_deer %>%
    dplyr::select(Species, SiteID, Distance, size, Date, Time, Behaviour) %>%
    dplyr::full_join(cams_curated2 %>%
                       dplyr::select(SiteID, Tkt, Effort), by="SiteID") %>%
    dplyr::mutate(object=1:nrow(.)) %>%
    dplyr::mutate(size = if_else(is.na(size),0L, size)) %>%
    dplyr::arrange(SiteID)

  ### Format data in a time format for availability
  summarised_count <- dcount
  summarised_count$Distance <- factor(summarised_count$Distance, levels = c("0 - 2.5", "2.5 - 5", "5 - 7.5", "7.5 - 10", "10+"))
  
  summarised_count <- summarised_count %>%
    mutate(distance = case_when(Distance == "0 - 2.5" ~ 1.25,
                                Distance == "2.5 - 5" ~ 3.75,
                                Distance == "5 - 7.5" ~ 6.25,
                                Distance == "7.5 - 10" ~ 8.75,
                                Distance == "10+" ~ 11.25)) %>%
    filter(size == 1)
  
  # site variables 
    site_vars <- cams_curated2 %>%
    left_join(dplyr::tbl(con, dbplyr::in_schema("deervic", "curated_site_data")) %>%
    dplyr::filter(SiteID %in% !!c(cams_curated2$SiteID, "47191")) %>%
    dplyr::collect() %>%
    dplyr::mutate(HerbaceousUnderstoryCover = scale(NNWHUCover + ENWHUCover),
                  WoodyUnderstoryCover = scale(NWUCover + EWUCover),
                  Trail = factor(Trail),
                  SiteID = dplyr::case_when(SiteID == "47191" & CameraID == "HO04101053" ~ "47191A",
                                            SiteID == "47191" & CameraID != "HO04101053" ~ "47191B",
                                            TRUE ~ SiteID))) %>% # native + exotic herbaceous cover
    dplyr::arrange(SiteID) %>%
    as.data.frame() %>%
      select(-Time, -Tkt, -Date, -Effort) 
    
    ds_data <- summarised_count %>% left_join(site_vars) %>%
      filter(!is.na(HerbaceousUnderstoryCover))
  
  mybreaks <- c(0,2.5,5,7.5,10,12.5)
  trunc.list <- list(left=0, right=12.5)
  conversion <- convert_units("metre", "metre", "metre")
  
  hn0 <- ds(ds_data, transect = "point", key="hn", adjustment = NULL,
          convert_units = conversion, cutpoints = mybreaks, truncation = trunc.list)
  
  hr0 <- ds(ds_data, transect = "point", key="hr", adjustment = NULL,
          convert_units = conversion, cutpoints = mybreaks, truncation = trunc.list)
  
  hn1 <- ds(ds_data, transect = "point", key="hn", adjustment = NULL, 
            formula = ~HerbaceousUnderstoryCover,
          convert_units = conversion, cutpoints = mybreaks, truncation = trunc.list)
  
  hr1 <- ds(ds_data, transect = "point", key="hr", adjustment = NULL, 
            formula = ~HerbaceousUnderstoryCover,
          convert_units = conversion, cutpoints = mybreaks, truncation = trunc.list)
  
  kableExtra::kbl(summarize_ds_models(hn0, hr0, hn1, hr1, output = "latex") %>% 
                    mutate(Model = stringr::str_extract(Model, pattern = '(?<=\\{)[^\\}]+')), 
                  digits=3,
             caption="Model selection table for deer detection", format = "html") %>%
    kable_styling(bootstrap_options = "striped")
Table 1: Model selection table for deer detection
Model Key function Formula \(\chi^2\) \(p\)-value \(\hat{P_a}\) se(\(\hat{P_a}\)) \(\Delta\)AIC
4 hr1 Hazard-rate ~HerbaceousUnderstoryCover 0 0.230 0.008 0.000
2 hr0 Hazard-rate ~1 0 0.228 0.008 3.759
1 hn0 Half-normal ~1 0 0.336 0.004 186.719
3 hn1 Half-normal ~HerbaceousUnderstoryCover 0 0.336 0.004 187.700

Model Options

Below we list the models we are implementing, all have a hazard detection function, but we compare four different abundance formulas as well as the inclusion of the EVC random effect.

Show code
#Models to run 
models_to_run <- data.frame(formula = c("ab_formula_1", 
                                        "ab_formula_2",
                                        "ab_formula_3",
                                        "ab_formula_4",
                                        "ab_formula_3", 
                                        "ab_formula_4",
                                        "ab_formula_1",
                                        "ab_formula_2"), 
                            keyfun = c(1,1,1,1,1,1,1,1), 
                            evc = c(F,F,F,F,T,T,T,T))

Full Models

We can fit models using cmdstanr. Here we fit models using a poisson distribution. The models we compare are:

Show code
model_fits <- list()

for(i in 1:nrow(models_to_run)) {

  if(models_to_run$evc[i]) {
    model_to_fit <- model_negbin_ms_evc
  } else {
    model_to_fit <- model_negbin_ms
  }
  
  form_to_use <- get(models_to_run$formula[i])
  
  params_to_use <- c(TRUE, labels(terms(ab_formula_1)) %in% labels(terms(form_to_use)))
  
  data_to_use <- multispecies_data
  
  data_to_use$X_psi <- as.matrix(data_to_use$X_psi[,params_to_use])
  data_to_use$X_pred_psi <- as.matrix(data_to_use$X_pred_psi[,params_to_use])
  data_to_use$m_psi <- sum(params_to_use)
  
  data_to_use$keyfun <- models_to_run$keyfun[i]

  model_fits[[i]] <- model_to_fit$sample(data = data_to_use,
                                         chains = nc,
                                 parallel_chains = nc, 
                                 init = 0.1, 
                                 max_treedepth = 10, 
                                 refresh = 25, 
                                 show_messages = TRUE,
                                 show_exceptions = FALSE,
                                 save_warmup = FALSE,
                                 iter_sampling = ni, 
                                 iter_warmup = nw)

  model_fits[[i]]$save_object(paste0("outputs/models/fit_multispecies_", 
                                     i,".rds"))

  
}
Show code
prefix <- "outputs/models/fit_multispecies_"
model_fits <- list()
for(i in 1:nrow(models_to_run)) {
  file_to_read <- paste0(prefix, i, ".rds")
  
  model_fits[[i]] <- readRDS(file_to_read)
}

Model Evaluation

Our strategy for model evaluation is to determine the best performing model out of a limited set of options with different fixed-effects, random-effects and detection functions.

LOO (Leave-One-Out Cross-Validation).

We use leave-one-out cross-validation (LOO-CV); a Bayesian model evaluation process (Vehtari, Gelman, and Gabry 2017; Vehtari et al. 2020) to compare the models. Better performing models according to loo::loo() will have higher elpd values.

Show code
loos <- list()
for(i in 1:length(model_fits)) {
  loos[[i]] <- model_fits[[i]]$loo(variables = "log_lik2_species", cores = 8)
}

names(loos) <- paste("model", 1:length(model_fits))

loo_compare_table <- loo::loo_compare(x = loos)

# Plot loo table
  gt(loo_compare_table %>% 
       as.data.frame() %>% 
              tibble::rownames_to_column()) %>% 
    fmt_number(decimals = 2) %>%
    tab_style(locations = cells_row_groups(), style = cell_text(weight = "bold"))
elpd_diff se_diff elpd_loo se_elpd_loo p_loo se_p_loo looic se_looic
model 7 0.00 0.00 2,120.29 7.66 1.24 0.06 −4,240.59 15.32
model 8 −0.09 0.30 2,120.20 7.64 1.19 0.05 −4,240.41 15.29
model 5 −0.20 0.32 2,120.09 7.64 1.17 0.05 −4,240.19 15.27
model 6 −1.76 0.68 2,118.53 7.66 1.04 0.05 −4,237.06 15.32
model 1 −1.94 0.62 2,118.35 7.65 1.04 0.05 −4,236.70 15.29
model 2 −2.16 0.68 2,118.13 7.65 1.00 0.05 −4,236.26 15.30
model 3 −2.40 0.68 2,117.89 7.64 0.98 0.05 −4,235.79 15.28
model 4 −4.37 0.92 2,115.93 7.69 0.85 0.04 −4,231.86 15.38

Given that models are similar in there predictive performance at sampled sites we also evaluate the performance of predictions across unsampled sites. Below we see that more complex models (usually having a slighly better elpd), have poorer performance to unsampled areas and lead to unrealistic predictions.

Show code
Nhat_CV <- list()
for(i in 1:length(model_fits)) {
  Nhat_CV[[i]] <- model_fits[[i]]$summary("Nhat", CV = ~ sd(.x)/mean(.x)) %>%
    mutate(Model = i, Species = deer_species_all)
}

Nhat_CV_c <- bind_rows(Nhat_CV)

Nhat_CV_c %>%
  select(Model, Species, CV) %>%
  kbl("html", escape=FALSE, caption = "Coefficients of Variation for Total Abundance estimates for the models assessed in the analysis") %>%
  kable_styling() %>%
  kableExtra::column_spec(3, background = ifelse(Nhat_CV_c$CV > 0.5, "red", ifelse(Nhat_CV_c$CV > 0.2, "orange", "green")), color = "white") %>%
 collapse_rows(1) 
Table 2: Coefficients of Variation for Total Abundance estimates for the models assessed in the analysis
Model Species CV
1 Cervus unicolor 0.1886410
Dama dama 0.9422867
Cervus elaphus 3.4679094
Axis porcinus 0.6694140
2 Cervus unicolor 0.1595965
Dama dama 0.4671161
Cervus elaphus 1.1285972
Axis porcinus 0.5351523
3 Cervus unicolor 0.1589921
Dama dama 0.4140389
Cervus elaphus 0.8865843
Axis porcinus 0.5135488
4 Cervus unicolor 0.1516583
Dama dama 0.2145688
Cervus elaphus 0.6266303
Axis porcinus 0.3045017
5 Cervus unicolor 0.2132886
Dama dama 0.5446794
Cervus elaphus 3.3479264
Axis porcinus 0.7870094
6 Cervus unicolor 0.2062165
Dama dama 0.2894896
Cervus elaphus 1.5069138
Axis porcinus 0.4926303
7 Cervus unicolor 0.2586611
Dama dama 0.5460739
Cervus elaphus 4.9932521
Axis porcinus 0.9147719
8 Cervus unicolor 0.2183291
Dama dama 0.4281970
Cervus elaphus 4.1957471
Axis porcinus 1.0076050

Selecting the ‘top’ model

Based on the findings of loo and comparisons of other checks (posterior predictive checks, variation in posterior draws, plausibility of predictions), we see that all model perform similarly well using loo. However, several models produce predictions with lower levels of variation in the posterior predictions. Based on this, we use model 3 hereafter for analysis. This model has fixed effects of scale(sqrt(BareSoil)), scale(sqrt(Nitrogen)), scale(log(PastureDistance)), scale(BIO15) and scale(sqrt(ForestEdge)). It also has a random effect of bioregion for each species and an over-dispersion parameter in the camera counts based on bioregion. This model does not have as many fixed effects as models 1 and 2, and also does not have a second random effect of EVC alongside bioregion (model 6). Unfortunately, the additional random effect of EVC produces much more variable and over-dispersed results for Red Deer, so we opt to use model 3 instead.

Show code
# assign index for top model
top <- 3

Posterior predictive checks

Posterior predictive checks allow us to compare the observed data to the model-generated data. For each species we undertake posterior predictive checks for summary statistics relating to the number of deer seen on the cameras at each site. Ideally a well-fit model is able to make predictions that match the observed data. Here counts are the number of snapshot moments with deer. The summary statistics we use for the posterior predictive checks are:

  1. Maximum counts of deer seen at a site on a camera
  2. Mean number of deer seen at a site on a camera
  3. Standard deviation of the counts of deer on the cameras
  4. Average (mean) counts at sites in a scatter plot
  5. Proportion of sites with zero counts (camera or transect)
  6. Observed vs expected proportion of counts

Sambar Deer

Based on the plot below, we see that the PPCs for Sambar perform sufficiently well.

Show code
q95 <- function(x) quantile(x, 0.95, na.rm = T)
# q25 <- function(x) quantile(x, 0.25, na.rm = T)
q75 <- function(x) quantile(x, 0.75, na.rm = T)
sd_90 <- function (x, na.rm = FALSE) {
  quants <- quantile(x, c(0.05, 0.95), na.rm = T)
  x <- x[x < quants[2] & x > quants[1]]
  sqrt(var(if (is.vector(x) || is.factor(x)) x else as.double(x),
    na.rm = na.rm))
}

bayes_theme <- c(delwp_cols_seq$Teal[6:8], delwp_cols_shades$Navy[1:3])

bayesplot::color_scheme_set(bayes_theme)

funs <- c(max, mean, sd, ppc_scatter_avg, prop_zero, ppc_dens_overlay)
titles <- c("Max", "Mean", "SD", "Average Site Counts", "Proportion Zeros", "Density of Counts")

ppc_sambar <- list()

# funs <- c(prop_zero, mean, q90, sd)

for(i in 1:length(funs)) {
ppc_sambar[[i]] <- posterior_checks_multispecies(model = model_fits[[top]], 
                 model_data = multispecies_data, species_index = 1,
                 stat = funs[[i]], integrated = T, only_det = F,
                 title = titles[i])
}

cowplot::plot_grid(plotlist = ppc_sambar, labels = "AUTO", ncol = 2)
Posterior predictive checks for Sambar Deer

Figure 1: Posterior predictive checks for Sambar Deer

Fallow Deer

The PPC’s below show that fallow has more dispersion in camera counts than Sambar. While the observed counts fall within the intervals of our model, there is higher uncertainty. Some of these findings appear to be driven by a site with very high counts (over 2000 snapshot moments of deer individuals), this is more than twice the second highest.

Show code
ppc_fallow <- list()

for(i in 1:length(funs)) {
ppc_fallow[[i]] <- posterior_checks_multispecies(model = model_fits[[top]], 
                 model_data = multispecies_data, species_index = 2,
                 stat = funs[[i]], integrated = T, only_det = F,
                 title = titles[i])
}

cowplot::plot_grid(plotlist = ppc_fallow, labels = "AUTO", ncol = 2)
Posterior predictive checks for Fallow Deer

Figure 2: Posterior predictive checks for Fallow Deer

Red Deer

Red Deer have PPC’s that show a high congruence between predicted and expected counts.

Show code
ppc_red <- list()

for(i in 1:length(funs)) {
ppc_red[[i]] <- posterior_checks_multispecies(model = model_fits[[top]], 
                 model_data = multispecies_data, species_index = 3,
                 stat = funs[[i]], integrated = T, 
                 title = titles[i])
}

cowplot::plot_grid(plotlist = ppc_red, labels = "AUTO", ncol = 2)
Posterior predictive checks for Red Deer

Figure 3: Posterior predictive checks for Red Deer

Hog Deer

Hog Deer have PPC’s that show our model performs sufficiently well in explaining the patterns of observed counts.

Show code
ppc_hog <- list()

for(i in 1:length(funs)) {
ppc_hog[[i]] <- posterior_checks_multispecies(model = model_fits[[top]], 
                 model_data = multispecies_data, species_index = 4,
                 stat = funs[[i]], integrated = T,
                 title = titles[i])
}

cowplot::plot_grid(plotlist = ppc_hog, labels = "AUTO", ncol = 2)
Posterior predictive checks for Hog Deer

Figure 4: Posterior predictive checks for Hog Deer

Model Convergence

We observed no divergences or any STAN sampling warnings/issues for our models. To visualise the convergence of the model we can observe the mixing of chains for key parameters below:

Show code
bayesplot::color_scheme_set(unname(delwp_cols[1:6]))
convergence_params <- model_fits[[top]]$draws(c("beta_det", 
                                                "beta_psi",
                                                "bioregion_sd", 
                                                "beta_trans_det", 
                                                "odRN", 
                                                "od_mu"))

mcmc_trace(convergence_params, facet_args = list(ncol = 4)) + 
  theme(panel.spacing = unit(0.5, "lines"))
Chain mixing of several key parameters relating to detection, abundance and dispersion

Figure 5: Chain mixing of several key parameters relating to detection, abundance and dispersion

Model Predictions

Within the STAN model we generate predictions for sampled and un-sampled locations. This provides us with site-level abundance estimates as well as estimates across all (un-sampled) public forest.

Site-based Predictions

Site-based Detections Map

Below we display a detection/non-detection map for all sites.

Show code
vic_regions <- vicmap_query("open-data-platform:delwp_region") %>%
  collect() %>%
  st_transform(3111) %>%
  st_simplify(dTolerance = 500)

delwp_pal <- colorRampPalette(c("#B9C600",
                                delwp_cols[["Teal"]], 
                                delwp_cols[["Navy"]]))

site_detection_plot <- function(data, regions, cams_curated, species_index) {
  
  cam_data <- cams_curated %>% dplyr::select(SiteID)
  
  cam_data$Detected <- factor(data$any_seen[species_index, ])
  cam_data$cam_seen <- factor(as.integer(as.logical(rowSums(data$n_obs[,,species_index]))))
  
  if(species_index == 4) {
    regions <- regions %>% filter(delwp_region == "GIPPSLAND")
    cam_data <- cam_data %>% st_filter(regions %>% st_transform(4283))
  }
  
   lvls <- length(unique(cam_data$Detected, na.rm = T))
  
plot <- ggplot2::ggplot(data = cam_data) +
  ggplot2::geom_sf(data = regions, alpha = 0.75, fill = "grey80") +
  ggplot2::geom_sf(aes(fill = Detected), shape = 21, size = 2, alpha = 0.75) +
  ggplot2::scale_fill_manual(values = delwp_pal(lvls), 
                             na.value = "transparent", na.translate = F,
                             name = "Detected", guide = guide_legend(override.aes = list(size = 6))) +
  # ggtitle(species_names[species_index])+
  theme_bw() +
  theme(legend.text = element_text(size = 18), legend.key.size = unit(1, "cm"),
        title = element_text(size = 20))
  
return(plot)
  
}

detection_plot <- cowplot::plot_grid(site_detection_plot(multispecies_data, regions = vic_regions, cams_curated = cams_curated, species_index = 1),
site_detection_plot(multispecies_data, regions = vic_regions, cams_curated = cams_curated, species_index = 2),
site_detection_plot(multispecies_data, regions = vic_regions, cams_curated = cams_curated, species_index = 3),
site_detection_plot(multispecies_data, regions = vic_regions, cams_curated = cams_curated, species_index = 4), ncol = 2, 
labels = c("A) Sambar Deer", "B) Fallow Deer", "C) Red Deer", "D) Hog Deer"), greedy = TRUE)

ggsave(plot = detection_plot, filename = "outputs/plots/all_detection_plot.png", width = 4800, height = 3200, units = "px", bg = "white")

detection_plot
Detections of deer (camera or transects) from across Victoria. Detections of Hog Deer are restricted to Gippsland.

Figure 6: Detections of deer (camera or transects) from across Victoria. Detections of Hog Deer are restricted to Gippsland.

Site-based Prediction Map

Below we plot visualisations of point-estimates for the various deer species surveyed for in this study.

Show code
site_preds <- function(model, cams_curated, species_index) {
  
rn_dens <- model$summary("N_site", prob_occ = ~ mean(. > 0), trimmed_mean = ~mean(., trim = 0.025), CV = ~sd(.)/mean(.), posterior::default_summary_measures(), posterior::default_convergence_measures())
  
which_sp <- which(stringr::str_detect(string = rn_dens$variable,
                                         pattern = paste0("N_site\\[",
                                                          species_index)))
  

  density_at_sites_rn <- cbind(cams_curated, rn_dens[which_sp, ])
  
  return(density_at_sites_rn)
}

site_density_plot <- function(densities, regions, species) {
  densities$density <- cut(densities$mean,
                                   breaks = c(0, 0.25, 0.5, 1, 3, 5, 10, max(densities$mean)),
                                   labels = c("< 0.25",
                                                    "0.25 - 0.5", 
                                                    "0.5 - 1", 
                                                    "1 - 3", 
                                                    "3 - 5",
                                                    "5 - 10", 
                                                    "10 +"), include.lowest = T, right = T)
  
  if(species == "Hog") {
    regions <- regions %>% filter(delwp_region == "GIPPSLAND")
    densities <- densities %>% st_filter(regions %>% st_transform(4283))
  }
  
   lvls <- length(unique(densities$density, na.rm = T))
  
plot <- ggplot2::ggplot(data = densities) +
  ggplot2::geom_sf(data = regions, alpha = 0.75, fill = "grey80") +
  ggplot2::geom_sf(aes(fill = density, alpha = mean), shape = 21, size = 3) +
  ggplot2::scale_fill_manual(values = delwp_pal(lvls), 
                             na.value = "transparent", na.translate = F,
                             name = "", guide = guide_legend(override.aes = list(size = 6))) +
  scale_alpha_continuous(range = c(0.5,1), guide = "none") +
  labs(title = paste0('Abundance of ',species ,' Deer'), fill = bquote('Deer per'~km^2)) +
  theme_bw() +
  theme(legend.text = element_text(size = 18), legend.key.size = unit(1, "cm"),
        title = element_text(size = 20))
  
return(plot)
  
}

sambar_preds <- site_preds(model_fits[[top]], cams_curated = cams_curated, species_index = 1)
fallow_preds <- site_preds(model_fits[[top]], cams_curated = cams_curated, species_index = 2)
red_preds <- site_preds(model_fits[[top]], cams_curated = cams_curated, species_index = 3)
hog_preds <- site_preds(model_fits[[top]], cams_curated = cams_curated, species_index = 4)

site_deer_predictions <- list(Sambar = sambar_preds, 
                         Fallow = fallow_preds, 
                         Red = red_preds, 
                         Hog = hog_preds)

saveRDS(site_deer_predictions, "outputs/site_deer_predictions.rds")

cowplot::plot_grid(site_density_plot(sambar_preds, vic_regions, species = "Sambar"),
site_density_plot(fallow_preds, vic_regions, species = "Fallow"),
site_density_plot(red_preds, vic_regions, species = "Red"), 
site_density_plot(hog_preds, vic_regions, species = "Hog"), ncol = 1)
Site-level density estimates for all sites sampled as part of statewide and hog deer surveys. Point-density estimates are from trimmed means of the posterior draws (trimming of the top and bottom 1% of draws)

Figure 7: Site-level density estimates for all sites sampled as part of statewide and hog deer surveys. Point-density estimates are from trimmed means of the posterior draws (trimming of the top and bottom 1% of draws)

Average density estimates varied across sites for all four species. Sambar deer were the most commonly detected species on the camera traps (104/317 sites), and density estimates at sites ranged from 0 to 28.31 deer/km2, for fallow, fewer site detections (30), but more variability in site density gave a range between 0 and 23.25 deer/km2. Red deer were only detected on 12 sites with densities ranging from 0 to 6.81. Finally, Hog Deer, which were only detected on 22 cameras in South Gippsland, had site densities ranging from 0 to 16.9 deer/km2.

Site-density Summaries

As a sanity-check we compare the average modelled densities at sites with (i) no evidence of deer, (ii) evidence of deer present on transects, and (iii) evidence of deer present on cameras. We expect that average densities are generally higher at sites that detected some form of deer than sites that did not detect any sign of deer. Additionally, we would also expect average densities to be generally higher at sites that had detections on cameras than those with only detections from transects. The table below shows these expectations to be correct.

Show code
density_summary_table <- function(preds, model_data, species, species_index) {
  
  cam_seen <- as.integer(as.logical(rowSums(model_data$n_obs[,,species_index])))
  
  transect_seen <- model_data$transects %>%
    filter(Species == !!deer_species_all[species_index]) %>%
    mutate(NoCam = case_when(Survey == "Camera" ~ 0L, 
                             TRUE ~ 1)) %>%
    group_by(SiteID) %>%
    summarise(TransectDet = max(NoCam*Presence)) %>%
    pull(TransectDet)
  
  preds_sum <- preds %>% 
    st_drop_geometry() %>%
    mutate(`Species` = species,
           `CameraDetection` = cam_seen, 
           `AnyDetection` = model_data$any_seen[species_index, ], 
           `OnlyTransect`  = transect_seen,
           `Detection` = case_when(CameraDetection == 0 & AnyDetection == 0 ~ "Not seen", 
                                   CameraDetection == 0 & AnyDetection == 1 ~ "Only detected on transects", 
                                   CameraDetection == 1 & OnlyTransect == 0 ~ "Only seen on cameras", 
                                   CameraDetection == 1 & OnlyTransect == 1 ~ "Seen on both camera and transect")) %>%
    group_by(`Species`, `Detection`) %>%
    summarise(`Number of Sites` = n(),
              `Average Density` = mean(mean)) %>%
    ungroup()
  
  return(preds_sum)
}

density_summary <- bind_rows(
  density_summary_table(sambar_preds, multispecies_data, 
                        species = "Sambar", species_index = 1),
  density_summary_table(fallow_preds, multispecies_data, 
                        species = "Fallow", species_index = 2),
  density_summary_table(red_preds, multispecies_data, 
                        species = "Red", species_index = 3),
  density_summary_table(hog_preds, multispecies_data, 
                        species = "Hog", species_index = 4))

density_summary %>%
  kbl(format = "html", digits = 2, caption = "Average (mean) density estimates at the various sites groups of sites") %>%
  kable_styling("striped") %>%
  collapse_rows(1)
Table 3: Average (mean) density estimates at the various sites groups of sites
Species Detection Number of Sites Average Density
Sambar Not seen 182 0.81
Only detected on transects 31 2.46
Only seen on cameras 50 3.34
Seen on both camera and transect 54 4.16
Fallow Not seen 259 0.52
Only detected on transects 28 0.97
Only seen on cameras 20 1.37
Seen on both camera and transect 10 1.89
Red Not seen 304 0.05
Only detected on transects 1 1.23
Only seen on cameras 5 1.47
Seen on both camera and transect 7 2.01
Hog Not seen 293 0.27
Only detected on transects 2 2.03
Only seen on cameras 22 3.10

Regional and Statewide Abundance

Within our model we calculate abundance/density for deer in each of the 7.457^{4} km2 of public land. Based on these spatial predictions we can estimate abundance at a regional level (6 DEECA regions) and across the whole state.

Statewide Maps

Using the model predictions (“pred”) for all suitable public land we generate a raster (1km2 resolution). We save the average spatial estimates under outputs/rasters (tif files) and also provide binned plots (png files) in outputs/plots.

Show code
pred_raster_full <- terra::rast(prediction_raster)

pred_raster <- terra::app(pred_raster_full[[stringr::str_subset(
  stringr::str_remove_all(labels(terms(ab_formula_1)),
                          "scale[(]|[)]|log[(]|sqrt[(]"), 
  pattern = "[*]", negate = T)]], mean)

mean_raster <- list()
occ_raster <- list()
mean_raster_discrete <- list()

# gp_preds_draws_all <- model_fits[[top]]$draws("pred", format = "matrix")
gp_preds_draws_all <- model_fits[[top]]$draws("pred", format = "matrix")
# gp_preds_summary_all <- model_fits[[top]]$summary("pred", prob_occ = ~ mean(. > 0))


for(i in 1:length(deer_species_all)) {

which_sp <- which(stringr::str_detect(string = colnames(gp_preds_draws_all),
                                         pattern = paste0("pred\\[", i)))

gp_preds_draws_sp <- gp_preds_draws_all[,which_sp]

terra::values(pred_raster)[!is.na(terra::values(pred_raster))] <- apply(gp_preds_draws_sp, 
                                                                        2, 
                                                                        mean, 
                                                                        na.rm = T) #filter out highest draw

mean_raster[[i]] <- pred_raster
mean_raster_discrete[[i]] <- mean_raster[[i]]
max_pred <- max(values(mean_raster[[i]]), na.rm = T)

values(mean_raster_discrete[[i]]) <- cut(values(mean_raster_discrete[[i]]) , 
                                         breaks = c(0,0.25,0.5,1,3,5,10,max_pred), 
                                         include.lowest = T, right = T,
                                         labels = c("< 0.25",
                                                    "0.25 - 0.5", 
                                                    "0.5 - 1", 
                                                    "1 - 3", 
                                                    "3 - 5",
                                                    "5 - 10", 
                                                    "10 +"))

occ_raster[[i]] <- pred_raster
terra::values(occ_raster[[i]])[!is.na(terra::values(occ_raster[[i]]))] <- apply(gp_preds_draws_sp, 2, function(x) sum(x > 0)/length(x))

}

# combine mean rasters together
combined_raster <- rast(mean_raster)
names(combined_raster) <- c("Average Sambar Deer Density (per km2)", 
                            "Average Fallow Deer Density (per km2)", 
                            "Average Red Deer Density (per km2)", 
                            "Average Hog Deer Density (per km2)")

combined_occupancy_raster <- rast(occ_raster)
names(combined_occupancy_raster) <- c("Average Sambar Deer Occupancy Probability", 
                            "Average Fallow Deer Occupancy Probability", 
                            "Average Red Deer Occupancy Probability", 
                            "Average Hog Deer Occupancy Probability")

writeRaster(combined_occupancy_raster, "outputs/rasters/combined_occupancy_raster.tif", overwrite = T)
writeRaster(combined_raster, "outputs/rasters/combined_deer_average_density.tif", overwrite = T)
Show code
# reg <- VicmapR::vicmap_query("open-data-platform:delwp_region") %>%
#     VicmapR::collect() %>%
#   sf::st_transform(3111) %>%
#   sf::st_simplify(dTolerance = 250) 
state <- VicmapR::vicmap_query("open-data-platform:vmlite_victoria_polygon_su5") %>%
  filter(state != "NSW" & state != "SA" & feature_type_code != "sea") %>%
  VicmapR::collect() %>%
  sf::st_transform(3111)

gippsland <- vic_regions %>% filter(delwp_region == "GIPPSLAND")

plot_abundance <- function(raster, state, species, crop = NULL) {
  
  if(!is.null(crop)) {
    raster <- terra::crop(raster, vect(crop), mask = T)
    state <- crop
  }
  
  lvls <- length(unique(values(raster, na.rm = T)))
  
  plot <- ggplot2::ggplot() +
    ggplot2::geom_sf(data = state, alpha = 1, linewidth = 0.5, fill = "grey90") + 
    tidyterra::geom_spatraster(data = raster, na.rm = T) + 
    # tidyterra::scale_fill_terrain_d(na.translate = FALSE) + 
    ggplot2::scale_fill_manual(values = delwp_pal(lvls), na.value = "transparent", na.translate = F) +
    ggplot2::labs(fill = bquote('Deer per'~km^2)) + 
    ggspatial::annotation_scale() +
    delwp_theme()
  
  return(plot)
  
} 

sambar_abundance_plot <- plot_abundance(mean_raster_discrete[[1]], 
                                        state = state, 
                                        species = "Sambar Deer")
fallow_abundance_plot <- plot_abundance(mean_raster_discrete[[2]], 
                                        state = state, 
                                        species = "Fallow Deer")
red_abundance_plot <- plot_abundance(mean_raster_discrete[[3]], 
                                     state = state, 
                                     species = "Red Deer")
hog_abundance_plot <- plot_abundance(mean_raster_discrete[[4]], 
                                     state = state, crop = gippsland,
                                     species = "Hog Deer")

ggsave(plot = sambar_abundance_plot, filename = "outputs/plots/sambar_abundance_plot.png", 
       width = 2400, height = 1600, units = "px")
ggsave(plot = fallow_abundance_plot, filename = "outputs/plots/fallow_abundance_plot.png", 
       width = 2400, height = 1600, units = "px")
ggsave(plot = red_abundance_plot, filename = "outputs/plots/red_abundance_plot.png", 
       width = 2400, height = 1600, units = "px")
ggsave(plot = hog_abundance_plot, filename = "outputs/plots/hog_abundance_plot.png", 
       width = 2400, height = 1600, units = "px")

ggsave(plot = cowplot::plot_grid(sambar_abundance_plot, fallow_abundance_plot, red_abundance_plot, hog_abundance_plot, ncol = 2, labels = c("A) Sambar Deer", "B) Fallow Deer", "C) Red Deer", "D) Hog Deer"), greedy = F), filename = "outputs/plots/all_abundance_plot.png", width = 4800, height = 3200, units = "px", bg = "white")

cowplot::plot_grid(sambar_abundance_plot, fallow_abundance_plot, red_abundance_plot, hog_abundance_plot, ncol = 1, labels = c("A) Sambar Deer", "B) Fallow Deer", "C) Red Deer", "D) Hog Deer"))
Abundance of Sambar, Fallow, Red and Hog Deer across Victoria, dark-grey area reflects area not included in predictions (i.e. not public land)

Figure 8: Abundance of Sambar, Fallow, Red and Hog Deer across Victoria, dark-grey area reflects area not included in predictions (i.e. not public land)

Regional Abundance Estimates

For each DEECA region we provide species-level estimates of abundance with 90 % confidence interval. We also calculate the model-based average density within each region based on the abundance and the total area of public land within each region.

Show code
Nhat_sum <- model_fits[[top]]$summary("Nhat", trimmed_mean = ~mean(., trim = 0.025), CV = ~sd(.)/mean(.), posterior::default_summary_measures(), posterior::default_convergence_measures())

Nhat_all <- model_fits[[top]]$summary("Nhat_sum", trimmed_mean = ~mean(., trim = 0.025), CV = ~sd(.)/mean(.), posterior::default_summary_measures(), posterior::default_convergence_measures())

format_nhat <- function(sp, data = Nhat_sum) {
  fd <- data %>%
    filter(variable == !!paste0("Nhat[", sp, "]"))
  
  paste0("$\\hat N$ = ",
         formatC(fd$trimmed_mean, digits = 0, big.mark = ",", format = "f"), 
        " [90% CI: ", 
        formatC(fd$q5, digits = 0, big.mark = ",", format = "f"), 
        " – ",
        formatC(fd$q95, digits = 0, big.mark = ",", format = "f"),
        "]")
}

reg <- VicmapR::vicmap_query("open-data-platform:delwp_region") %>%
    VicmapR::collect() %>%
    dplyr::group_by(delwp_region) %>%
    dplyr::summarise(geometry = sf::st_combine(geometry)) %>%
    dplyr::ungroup() %>%
    dplyr::mutate(delwp_region_fact = as.integer(factor(delwp_region))) 

abundance_table <- function(model, regions, pred_area, caption, species_index, return_data = F) {
  
  reg <- model$summary("Nhat_reg", prob_occ = ~ mean(. > 0), trimmed_mean = ~mean(., trim = 0.025), posterior::default_summary_measures(), posterior::default_convergence_measures())
  
  which_sp <- which(stringr::str_detect(string = reg$variable,
                                         pattern = paste0("Nhat_reg\\[", species_index)))
  
  regional_abundance <- reg[which_sp,] %>%
    dplyr::bind_rows(model$summary("Nhat", prob_occ = ~ mean(. > 0), trimmed_mean = ~mean(., trim = 0.025), posterior::default_summary_measures(), posterior::default_convergence_measures())[species_index,]) %>% 
    dplyr::mutate(variable = c(regions$delwp_region, "TOTAL"), 
                  `Area km2` = round(c(pred_area, sum(pred_area))), 
                  `Density (mean)` = round(trimmed_mean/`Area km2`, 2), 
                  `Density (5 %)` = round(q5/`Area km2`, 2),
                  `Density (95 %)` = round(q95/`Area km2`, 2)) %>% 
    dplyr::transmute(Region = variable, 
                  `Trimmed Mean` = round(trimmed_mean), 
                  Median = round(median), 
                  SD = round(sd), 
                  MAD = round(mad),
                  `5 %` = round(q5), 
                  `95 %` = round(q95), 
                  `Area km2`,
                  `Density (mean)`,
                  `Density (5 %)`,
                  `Density (95 %)`)
  
  if(return_data) {
    
 if(species_index == 4) {
    regional_abundance <- regional_abundance %>% 
      filter(Region == "GIPPSLAND")
  }
    
    return(regional_abundance %>% 
             mutate(Species = species_names[species_index], across(2:8, ~ formatC(.x, digits = 0, format = "f", big.mark = ",") )) %>%
             select(Species, everything()))
  }
  
  kableExtra::kbl(regional_abundance, digits = 2, format = "html", caption = caption, format.args = list(big.mark = ",")) %>%
    kableExtra::kable_styling("striped") %>%
    kableExtra::column_spec(1, bold = TRUE) %>%
    kableExtra::row_spec(6, hline_after = T) %>%
    kableExtra::row_spec(7, background = "#A5A1B5", color = "black", bold = T, hline_after = T)
}

pred_a_df <- data.frame(pred_reg = multispecies_data$pred_reg, 
                        area = multispecies_data$prop_pred) %>%
  group_by(pred_reg) %>%
  summarise(area = sum(area))

# Save Abundance Table 
all_ab_tab <- bind_rows(abundance_table(model_fits[[top]], regions = reg, pred_area = pred_a_df$area, species_index = 1, return_data = T), 
                        abundance_table(model_fits[[top]], regions = reg, pred_area = pred_a_df$area, species_index = 2, return_data = T),
                        abundance_table(model_fits[[top]], regions = reg, pred_area = pred_a_df$area, species_index = 3, return_data = T), 
                        abundance_table(model_fits[[top]], regions = reg, pred_area = pred_a_df$area, species_index = 4, return_data = T))

write.csv(all_ab_tab, "outputs/tables/regional_abundance_estimates.csv")

abundance_table(model_fits[[top]], regions = reg, pred_area = pred_a_df$area, caption = "Regional estimates of Sambar Deer Density", species_index = 1)
Table 4: Regional estimates of Sambar Deer Density
Region Trimmed Mean Median SD MAD 5 % 95 % Area km2 Density (mean) Density (5 %) Density (95 %)
BARWON SOUTH WEST 139 107 277 72 35 430 4,745 0.03 0.01 0.09
GIPPSLAND 54,895 54,109 8,898 8,675 42,618 70,948 25,211 2.18 1.69 2.81
GRAMPIANS 1,619 1,461 2,032 648 716 3,266 9,896 0.16 0.07 0.33
HUME 59,732 58,676 12,360 11,614 43,334 82,188 16,890 3.54 2.57 4.87
LODDON MALLEE 1,094 940 1,761 429 460 2,589 15,597 0.07 0.03 0.17
PORT PHILLIP 5,200 5,041 1,334 1,258 3,367 7,694 2,231 2.33 1.51 3.45
TOTAL 123,061 121,200 19,657 18,946 96,200 157,638 74,570 1.65 1.29 2.11
Show code
abundance_table(model_fits[[top]], regions = reg, pred_area = pred_a_df$area, caption = "Regional estimates of Fallow Deer Density", species_index = 2)
Table 4: Regional estimates of Fallow Deer Density
Region Trimmed Mean Median SD MAD 5 % 95 % Area km2 Density (mean) Density (5 %) Density (95 %)
BARWON SOUTH WEST 10,514 8,146 14,091 5,142 3,049 31,646 4,745 2.22 0.64 6.67
GIPPSLAND 12,305 11,592 5,001 3,711 6,914 21,794 25,211 0.49 0.27 0.86
GRAMPIANS 4,350 3,865 3,203 1,529 2,106 9,264 9,896 0.44 0.21 0.94
HUME 16,021 15,179 5,519 4,475 9,362 26,721 16,890 0.95 0.55 1.58
LODDON MALLEE 3,296 2,922 2,198 1,380 1,384 7,277 15,597 0.21 0.09 0.47
PORT PHILLIP 1,501 1,378 767 533 761 2,879 2,231 0.67 0.34 1.29
TOTAL 48,932 45,609 20,900 13,463 29,888 85,063 74,570 0.66 0.40 1.14
Show code
abundance_table(model_fits[[top]], regions = reg, pred_area = pred_a_df$area, caption = "Regional estimates of Red Deer Density", species_index = 3)
Table 4: Regional estimates of Red Deer Density
Region Trimmed Mean Median SD MAD 5 % 95 % Area km2 Density (mean) Density (5 %) Density (95 %)
BARWON SOUTH WEST 5,645 4,063 8,018 2,684 1,559 19,399 4,745 1.19 0.33 4.09
GIPPSLAND 179 143 196 92 49 502 25,211 0.01 0.00 0.02
GRAMPIANS 5,134 4,093 5,709 2,228 1,705 14,695 9,896 0.52 0.17 1.48
HUME 1,126 971 739 530 420 2,638 16,890 0.07 0.02 0.16
LODDON MALLEE 103 24 1,524 30 1 794 15,597 0.01 0.00 0.05
PORT PHILLIP 31 24 83 18 5 103 2,231 0.01 0.00 0.05
TOTAL 12,672 10,135 12,275 4,967 4,719 35,465 74,570 0.17 0.06 0.48
Show code
abundance_table(model_fits[[top]], regions = reg, pred_area = pred_a_df$area, caption = "Regional estimates of Hog Deer Density", species_index = 4) 
Table 4: Regional estimates of Hog Deer Density
Region Trimmed Mean Median SD MAD 5 % 95 % Area km2 Density (mean) Density (5 %) Density (95 %)
BARWON SOUTH WEST 0 0 0 0 0 0 4,745 0.00 0.00 0.00
GIPPSLAND 4,243 3,884 2,263 1,529 2,121 8,464 25,211 0.17 0.08 0.34
GRAMPIANS 0 0 0 0 0 0 9,896 0.00 0.00 0.00
HUME 0 0 0 0 0 0 16,890 0.00 0.00 0.00
LODDON MALLEE 0 0 0 0 0 0 15,597 0.00 0.00 0.00
PORT PHILLIP 0 0 0 0 0 0 2,231 0.00 0.00 0.00
TOTAL 4,243 3,884 2,263 1,529 2,121 8,464 74,570 0.06 0.03 0.11

Deer were found to be widely distributed across Victoria. Across approximately 7.457^{4} km2 of public land we estimate total deer abundance of the four species investigated in this study (Sambar, Fallow, Red and Hog) at 191,153 (90 % CI: 146,732 - 255,490. Sambar deer were the most populous species across Victoria (\(\hat N\) = 123,061 [90% CI: 96,200 – 157,638]), followed by Fallow at \(\hat N\) = 48,932 [90% CI: 29,888 – 85,063], Red (\(\hat N\) = 12,672 [90% CI: 4,719 – 35,465]), and Hog (\(\hat N\) = 4,243 [90% CI: 2,121 – 8,464]).

Threshold Distributions

From our abundance estimates we can obtain estimates of the distribution of the deer species across Victoria.

Show code
inv.cloglog <- function(x) 1-exp(-x)

detections <- function(data, regions, cams_curated, species_index) {
  
  cam_data <- cams_curated %>% dplyr::select(SiteID)
  
  cam_data$Detected <- factor(data$any_seen[species_index, ])
  cam_data$cam_seen <- factor(as.integer(as.logical(rowSums(data$n_obs[,,species_index]))))
  
  # if(species_index == 4) {
  #   regions <- regions %>% filter(delwp_region == "GIPPSLAND")
  #   cam_data <- cam_data %>% st_filter(regions %>% st_transform(4283))
  # }
  
   return(cam_data)
  
}

joint_threshold <- list()

combined_threshold <- list()

for(i in 1:length(species_names)) {

dets <- detections(multispecies_data, 
                       vic_regions, 
                       cams_curated = cams_curated, 
                   species_index = i) %>% 
  st_make_valid() %>%
  st_transform(3111)

pointVals <- terra::extract(combined_raster[[i]],
               terra::vect(dets %>% st_buffer(1000)), fun = "mean", na.rm = T)[,2]

rocData <- dets[which(!is.na(pointVals)),] %>% 
  mutate(lambda = na.omit(pointVals), 
         occ = inv.cloglog(lambda))

roc <- roc(rocData$Detected, rocData$lambda)

ci_coords <- ci.coords(roc, "best", 
                       best.policy = "random", 
                       best.method = "closest.topleft",
                       progress = "none", 
                       conf.level = 0.9)

threshold_rast <- combined_raster[[i]]
threshold_rast_lci <- combined_raster[[i]]
threshold_rast_uci <- combined_raster[[i]]

threshold_rast[values(threshold_rast) < ci_coords$threshold[2]] <- 0
threshold_rast[values(threshold_rast) >= ci_coords$threshold[2]] <- 1

threshold_rast_lci[values(threshold_rast_lci) < ci_coords$threshold[3]] <- 0
threshold_rast_lci[values(threshold_rast_lci) >= ci_coords$threshold[3]] <- 1

threshold_rast_uci[values(threshold_rast_uci) < ci_coords$threshold[1]] <- 0
threshold_rast_uci[values(threshold_rast_uci) >= ci_coords$threshold[1]] <- 1

combined_threshold[[i]] <- c(threshold_rast, threshold_rast_lci, threshold_rast_uci) %>% `names<-`(c("Threshold Median", "Threshold LCI", "Threshold UCI"))

joint_threshold[[i]] <- app(combined_threshold[[i]], sum, na.rm = T)

d <- data.frame(id = 0:3, values = c(NA_character_, "UCI", "Median", "LCI"))

levels(joint_threshold[[i]]) <- d

}

joint_threshold <- rast(joint_threshold)

names(combined_threshold) <- species_names
names(joint_threshold) <- species_names

writeRaster(joint_threshold, "outputs/rasters/combined_threshold_raster.tif", overwrite = TRUE)


get_th_area <- function(stack) {
  med <- sum(values(stack[[1]], na.rm = T) * multispecies_data$prop_pred)
  LCI <- sum(values(stack[[2]], na.rm = T) * multispecies_data$prop_pred)
  UCI <- sum(values(stack[[3]], na.rm = T) * multispecies_data$prop_pred)
  paste0(formatC(med, format = "fg", big.mark = ","), "km^2^ [95% CI: ",
         formatC(LCI, format = "fg", big.mark = ","), ", ",
         formatC(UCI, format = "fg", big.mark = ","), "]")
}

plot_threshold <- function(raster, state, crop = NULL) {
  
  if(!is.null(crop)) {
    raster <- terra::crop(raster, vect(crop), mask = T)
    state <- crop
  }
  
  # lvls <- length(unique(values(raster, na.rm = T)))
  
  plot <- ggplot2::ggplot() +
    ggplot2::geom_sf(data = state, alpha = 1, linewidth = 0.5, fill = "grey98") + 
    tidyterra::geom_spatraster(data = raster, aes(alpha = after_stat(value), 
                                                  fill = after_stat(value)), na.rm = T) + 
    # tidyterra::scale_fill_terrain_d(na.translate = FALSE) + 
    ggplot2::scale_fill_manual(values = unname(rep(delwp_cols[2], 3)), 
                               na.value = "transparent", 
                               na.translate = F) +
    ggplot2::scale_alpha_manual(values = c(0.25, 0.65, 0.85), 
                                na.value = "transparent", 
                                na.translate = F) +
    ggplot2::labs(alpha = "Distribution \nRange", fill = "Distribution \nRange") +
    ggspatial::annotation_scale() +
    delwp_theme()
  
  return(plot)
  
} 

sambar_th <- plot_threshold(joint_threshold[[1]], 
               state = state)

fallow_th <- plot_threshold(joint_threshold[[2]], 
               state = state)

red_th <- plot_threshold(joint_threshold[[3]], 
               state = state)

hog_th <- plot_threshold(joint_threshold[[4]], 
               state = state, crop = gippsland)

threshold_combined_plot <- cowplot::plot_grid(sambar_th, fallow_th, red_th, hog_th, ncol = 2, labels = c("A) Sambar Deer", "B) Fallow Deer", "C) Red Deer", "D) Hog Deer"), greedy = F)

ggsave(plot = threshold_combined_plot, filename = "outputs/plots/all_threshold_plot.png", width = 4800, height = 3200, units = "px", bg = "white")

threshold_combined_plot
Threshold distributions of Sambar, Fallow, Red and Hog Deer across Victoria, white area reflects area not included in predictions (i.e. not public land), LCI and UCI confidence bounds refer to bootstrapped threshold confidence intervals (90%) for optimal sensitivity and specificity.

Figure 9: Threshold distributions of Sambar, Fallow, Red and Hog Deer across Victoria, white area reflects area not included in predictions (i.e. not public land), LCI and UCI confidence bounds refer to bootstrapped threshold confidence intervals (90%) for optimal sensitivity and specificity.

We obtained a threshold distribution based on the mean abundance values for each grid cell. Our range estimates suggest Sambar Deer have a range of 38,582km2 [95% CI: 33,542, 40,596]. For Fallow deer the range estimate has more variabiliy 24,901km2 [95% CI: 20,666, 30,423], this is similarly the case for Red Deer (9,403km2 [95% CI: 2,465, 10,454]). Lastly, Hog Deer have a more isolated range estimate (2,176km2 [95% CI: 1,542, 2,219]).

Comparison of range size to existing distributions

Using range size maps based on presence-only records and VBA elicitation we can generate estimated existing distributions for the deer species on public land:

Show code
existing_dist <- list()
existing_dist[["Sambar"]] <- st_read("data/existing_distributions/Sambar_Deer_Distributions_Final_20201203")
Reading layer `Sambar_Deer_Distributions_Final_20201203' from data source `/Users/justincally/Documents/Work/R Repositories/statewide-deer-analysis/data/existing_distributions/Sambar_Deer_Distributions_Final_20201203' 
  using driver `ESRI Shapefile'
Simple feature collection with 10 features and 16 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 142.8374 ymin: -38.86406 xmax: 149.9739 ymax: -35.80198
Geodetic CRS:  GDA94
Show code
existing_dist[["Fallow"]] <- st_read("data/existing_distributions/Fallow_Deer_Distributions_Final_20201203")
Reading layer `Fallow_Deer_Distributions_Final_20201203' from data source `/Users/justincally/Documents/Work/R Repositories/statewide-deer-analysis/data/existing_distributions/Fallow_Deer_Distributions_Final_20201203' 
  using driver `ESRI Shapefile'
Simple feature collection with 70 features and 1 field
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 140.9643 ymin: -38.9087 xmax: 149.0549 ymax: -33.98139
Geodetic CRS:  GDA94
Show code
existing_dist[["Red"]] <- st_read("data/existing_distributions/Red_Deer_Distributions_Final_20201203")
Reading layer `Red_Deer_Distributions_Final_20201203' from data source `/Users/justincally/Documents/Work/R Repositories/statewide-deer-analysis/data/existing_distributions/Red_Deer_Distributions_Final_20201203' 
  using driver `ESRI Shapefile'
Simple feature collection with 32 features and 1 field
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 140.9687 ymin: -38.64039 xmax: 149.0283 ymax: -36.04568
Geodetic CRS:  GDA94
Show code
existing_dist[["Hog"]] <- st_read("data/existing_distributions/Hog_Deer_Distribution_20160329")
Reading layer `Hog_Deer_Distribution_20160329' from data source 
  `/Users/justincally/Documents/Work/R Repositories/statewide-deer-analysis/data/existing_distributions/Hog_Deer_Distribution_20160329' 
  using driver `ESRI Shapefile'
Simple feature collection with 1 feature and 48 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 145.7291 ymin: -39.13665 xmax: 149.2845 ymax: -37.72223
Geodetic CRS:  GDA94
Show code
crown_land <- readRDS("data/crown_land.rds") 

# Get overlapping CA 
existing_dist_cl <- lapply(existing_dist, FUN = function(x) {
  d <- st_transform(x, 3111)
  st_intersection(d, crown_land) %>%
    st_union()
})

existing_ranges <- sapply(existing_dist_cl, st_area) %>% 
  units::set_units("m2") %>% 
  units::set_units("km2") %>%
  formatC(., digits = 0, format = "fg", big.mark = ",") %>% 
  as.character()

names(existing_ranges) <- species_names
existing_ranges
Sambar Deer Fallow Deer    Red Deer    Hog Deer 
   "43,779"    "10,259"     "2,630"     "1,147" 

Deer Abundance in High Conservation Areas

Below we show code for extracting average abundance and densities for areas listed under the National Parks Act 1975, and greater than 20km2.

Show code
parkres <- vicmap_query("parkres") %>%
  collect() %>%
  st_transform(3111) %>%
  st_filter(vic_regions)

  |                                                        
  |                                                  |   0%
  |                                                        
  |=========================                         |  50%
  |                                                        
  |==================================================| 100%
Show code
parkres_grouped <- parkres %>%
  group_by(name, name_short) %>%
  summarise(area = sum(areasqm)/1e6, 
            geometry = st_combine(geometry)) %>%
  filter(area > 20) %>%
  ungroup()

ab_e <- terra::extract(combined_raster, 
                terra::vect(parkres_grouped), 
                fun = "sum", method = "simple", na.rm = T, touches = TRUE) %>%
  `names<-`(c("ID", species_names))

# get abundance estimates
park_ab <- bind_cols(parkres_grouped, round(ab_e[,2:5],0)) %>%
  na.omit() %>%
  mutate(`All Deer` = `Sambar Deer` + `Fallow Deer` + `Red Deer` + `Hog Deer`)%>% 
  mutate(across(contains("Deer"), ~ round(.x/area,2), .names = "Density {.col}")) %>% 
  mutate(Sambar = paste0(formatC(`Sambar Deer`, format =  "fg", big.mark = ","), " (", `Density Sambar Deer`, ")"), 
         Fallow = paste0(formatC(`Fallow Deer`, format =  "fg", big.mark = ","), " (", `Density Fallow Deer`, ")"), 
         Red = paste0(formatC(`Red Deer`, format =  "fg", big.mark = ","), " (", `Density Red Deer`, ")"), 
         Hog = paste0(formatC(`Hog Deer`, format =  "fg", big.mark = ","), " (", `Density Hog Deer`, ")"), 
         All = paste0(formatC(`All Deer`, format =  "fg", big.mark = ","), " (", `Density All Deer`, ")")) %>%
  arrange(desc(`Density All Deer`))

park_averages <- park_ab %>% 
  mutate(name = paste0(name, " (", round(area, 0), "km2)")) %>%
  dplyr::select(Reserve = name, 
            Sambar, Fallow, Red, Hog, All) %>%
  st_drop_geometry() 

# write.csv(park_averages, "outputs/tables/park_averages.csv", row.names = F)

Model Covariates

Our model uses covariates to inform three separate sub-models:

The following section explores the effect of these covariates and sub-models on detection and abundance.

Detection Processes

At a site, there are two processes that may lead us to not record deer on a camera when in reality they occupy the site (whereby the site is the home range area surrounding the camera.) Firstly, deer may be available for detection and enter the sampling area in front of the camera. However, due to a function of their distance from the camera we fail to detect this individual. Secondly, they may occupy a site (known from transect signs) but never enter the camera field of view, in this case we estimate the various detection rates from the various methods. The detection rate of the camera in this method is regarded as the spatial availability of deer for camera trap sampling.

Distance-sampling

Below we show the average detection rate for a given group of deer in front of the camera (up to 12.5m).

Show code
det_rates <- model_fits[[top]]$summary("p") %>%
  mutate(var = stringr::str_extract(variable, "(?<=\\[).+?(?=\\])")) %>%
  separate(var, into = c("site", "Group Size"), sep = ",")

av_det_rates <- det_rates %>%
  group_by(`Group Size`) %>%
  summarise(mean = mean(mean)) %>%
  ungroup() %>% 
  transmute(`Group Size`, 
            `Average Site Detection Probability` = mean)


av_det_rates %>%
  bind_rows() %>%
  arrange(`Group Size`) %>%
  kbl("html", digits = 3, caption = "Average detection rates for different group sizes") %>%
  kable_styling("striped") %>%
  collapse_rows(1)
Table 5: Average detection rates for different group sizes
Group Size Average Site Detection Probability
1 0.297
2 0.379
3 0.439
4 0.485
5 0.523

We can generate detection function plots that show how detection rates fall-off over varying distances. Initial exploration showed that there may be some relationship between herbaceous understorey and detection probability for Sambar Deer. However, when pooled together we see that herbaceous understorey appears to have no effect on detection.

Show code
site_vars <- dplyr::tbl(con, dbplyr::in_schema("deervic", "curated_site_data")) %>%
    dplyr::filter(SiteID %in% !!c(cams_curated$SiteID, "47191")) %>%
    dplyr::collect() %>%
    dplyr::mutate(HerbaceousUnderstoryCover = NNWHUCover + ENWHUCover,
           SiteID = dplyr::case_when(SiteID == "47191" & CameraID == "HO04101053" ~ "47191A",
                              SiteID == "47191" & CameraID != "HO04101053" ~ "47191B",
                              TRUE ~ SiteID)) %>% # native + exotic herbaceous cover
    dplyr::arrange(SiteID) %>%
    as.data.frame()

  n_distance_bins <- 5
  delta <- 2.5
  midpts <- c(1.25, 3.75, 6.25, 8.75, 11.25)
  max_distance <- 12.5
  gs <- 1:5
  
  pa <- matrix(data = NA, nrow = max(gs), ncol = length(midpts))
  for(j in 1:max(gs)) {
    pa[j, ] <- sapply(midpts, function(x) {pr_closest(x-(0.5*delta),
                                                      x+(0.5*delta),
                                                      max_distance = max_distance,
                                                      n = j)})
  }


det_curve <- model_fits[[top]]$draws("DetCurve", format = "draws_matrix") %>%
  as.data.frame() %>%
  # head(250) %>%
  pivot_longer(cols = everything())

det_curve_wr <- det_curve %>%
  mutate(var = stringr::str_extract(name, "(?<=\\[).+?(?=\\])")) %>%
  separate(var, into = c("s", "Distance"), sep = ",") %>%
  mutate(Distance = as.numeric(Distance)-0.5)

det_vars_pred <- site_vars %>%
  mutate(s = as.character(1:nrow(.)),
         herbaceouslvl = cut(HerbaceousUnderstoryCover,
                             breaks = c(0, 25, 50, 75, 105), 
                             labels = c("0 - 25%", 
                                        "25 - 50%", 
                                        "50 - 75%", 
                                        "75 - 100%"),
                             include.lowest = T, right = FALSE))

det_curve_sum <- det_curve_wr %>%
  mutate(Distance = as.numeric(Distance)) %>%
  left_join(det_vars_pred) %>%
  group_by(herbaceouslvl, Distance) %>%
  summarise(median = quantile(value, 0.5),
            q5 = quantile(value, 0.05),
            q95 = quantile(value, 0.95))

scaleHN <- det_curve_sum %>% filter(Distance == 1.5) %>% pull(median) %>% mean() # Approx scaling

y_comb_list <- list()
for(i in 1:max(gs)) {
y_comb_list[[i]] <- colSums(multispecies_data$y[,,i]) %>% 
  as.data.frame() %>%
  `colnames<-`("Count") %>%
  mutate(Distance = midpts,
         Bin = 1:5,
         gs = i,
         Prop = Count/(max(Count)),
         CountS = Count/pa[Bin,i])

}

y_combined <- bind_rows(y_comb_list) %>% 
  group_by(Distance) %>%
  mutate(SumCount = sum(Count)) %>%
  ungroup() %>%
  mutate(propC = Count/SumCount, 
         PropSM = (CountS/max(CountS)) * propC)  %>%
  group_by(Distance) %>% 
  summarise(PropS = sum(PropSM) * scaleHN)

detection_plot_HN <- ggplot(aes(x = Distance), data = det_curve_sum) +
  geom_col(aes(y = PropS), fill = "grey70", colour = "grey20", width = 2.5, data = y_combined, alpha = 0.7) +
  # geom_error bar(aes(ymin = PropSq5, ymax = PropSq95),  data = y_combined) +
  # geom_line(aes(y = HNS)) +
  geom_ribbon(aes(ymin = q5, ymax = q95, fill = herbaceouslvl), alpha = 0.5) +
  geom_line(aes(y = median, colour = herbaceouslvl)) +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0), limits = c(0,1)) +
  scale_fill_manual(values = unname(delwp_cols[c(1,2,3,8)])) +
  scale_colour_manual(values = unname(delwp_cols[c(1,2,3,8)])) +
  ylab("Detection probability (p)") +
  labs(fill = "Herbaceous\nUnderstorey\nCover", colour = "Herbaceous\nUnderstorey\nCover") +
  theme_classic()

ggsave(plot = detection_plot_HN, filename = "outputs/plots/hazard_detection_plot.png", width = 2400, height = 1600, units = "px")

detection_plot_HN
Distance-sampling detection process for Deer (group size = 1-5). Herbaceous understorey cover in the area surrounding appears to have some impact on detection probability

Figure 10: Distance-sampling detection process for Deer (group size = 1-5). Herbaceous understorey cover in the area surrounding appears to have some impact on detection probability

Transect-surveys

For Sambar, Fallow and Red Deer we were able to estimate detection rates of the various methods used to detect deer (camera, footprints, pellets, rubbings and wallows).

Show code
all_draws <- model_fits[[top]]$draws("beta_trans_det", format = "df")

# det_marginal_effects <- list()
# det_plot <- list()

obs_vars <- unlist(stringr::str_remove_all(string = labels(terms(transect_formula)), pattern =  "log|scale|\\)|\\("))
obs_cols <- unlist(stringr::str_remove_all(string = colnames(multispecies_data$trans_pred_matrix), pattern =  "log|scale|\\)|\\("))
obs_logs <- unlist(stringr::str_detect(string = colnames(multispecies_data$trans_pred_matrix), pattern =  "log\\("))

obs_labs <- c("Survey Method")

fac <- c(TRUE, TRUE, TRUE, TRUE, TRUE)
fac2 <- c(TRUE)

det_obs <- multispecies_data$transects %>% 
  mutate(Survey = factor(Survey))

det_marginal_effects <- list()
det_plot <- list()

params_w_levs <- levels(det_obs$Survey)

for(i in 1:(length(obs_cols))) {
det_marginal_effects[[i]] <- marginal_effects_cmd(all_draws, 
                                              param = "beta_trans_det", 
                                              param_number = i, log = obs_logs[i],
                                     model_data = det_obs, 
                                     model_column = obs_vars[c(1,attr(multispecies_data$trans_pred_matrix, "assign")[-1])[i]], 
                                     transition = FALSE) %>%
  mutate(group = param, 
         param = params_w_levs[i], 
         factor = fac[i], 
         variable = as.numeric(variable))

if(fac[i]) {
  det_marginal_effects[[i]] <- det_marginal_effects[[i]]
}

}

marginal_prob <- function(x, pwr = 3) {
  xm <- 1-x
  return(1-(xm^pwr))
}

det_marginal_effects_bind <- bind_rows(det_marginal_effects) %>%
  rowwise() %>%
  mutate(value = case_when(param != "Camera" ~ marginal_prob(value), 
                           TRUE ~ value))

det_marginal_effects_split <- split(det_marginal_effects_bind, det_marginal_effects_bind$group)

for(i in 1:length(det_marginal_effects_split)) {
  
  if(fac2[i]) {
    plot_data <- det_marginal_effects_split[[i]] %>% 
      mutate(variable = param, 
             param = group)
  } else {
    plot_data <- det_marginal_effects_split[[i]]
  }

det_plot[[i]] <- plot_data

}
data_for_plot <- bind_rows(det_plot) %>% 
  mutate(species = "All (combined)") %>%
  mutate(Density = "1 deer/km^2^")

data_d3 <- data_for_plot %>% 
  mutate(value = 1 - (1 - value)^3, 
         Density = "3 deer/km^2^")

data_d5 <- data_for_plot %>% 
  mutate(value = 1 - (1 - value)^5, 
         Density = "5 deer/km^2^")

data_combined <- bind_rows(data_for_plot, data_d3, data_d5)

transect_det_plot <- marginal_effects_plot_cmd_all(data_combined, 
                              factor = TRUE,
                              ylab = "Detection probability") +
  xlab("Survey method") +
  scale_y_continuous(limits = c(0,1)) +
  theme(legend.position = "top") +
  scale_fill_manual(values = unname(delwp_cols[c(1,2,3)])) +
  # scale_x_discrete(expand = c(0, 0)) +
  theme_classic() +
  theme(legend.text = ggtext::element_markdown())

ggsave(plot = transect_det_plot, filename = "outputs/plots/transect_detection_plot.png", width = 2400, height = 1400, units = "px")

transect_det_plot
Unconditional detection probabilities for the various methods of survey, across various deer densities. Camera trap detection probability is based on average deployment length (53 days), while transects are based on 3 transect searches

Figure 11: Unconditional detection probabilities for the various methods of survey, across various deer densities. Camera trap detection probability is based on average deployment length (53 days), while transects are based on 3 transect searches

Abundance Processes

Bioregion

We used a spatially-derived random-effect of the bioregion of the sampled site. This random-effect allows us the ability to make predictions that include the variance associated with this random effect. This random-effect also minimises predictions in areas without deer detection’s (e.g. the mallee).

Show code
bioregion_sd <- model_fits[[top]]$summary("bioregion_sd")

bayesplot::color_scheme_set(bayes_theme)
bioregion_contribution <- function(model, data, species, species_index) {
  eps_bioregion <- model$draws("eps_bioregion", format = "matrix")
  
  which_sp <- which(stringr::str_detect(string = colnames(eps_bioregion),
                                         pattern = paste0("eps_bioregion\\[", i)))
  
  eps_bioregion <- eps_bioregion[,which_sp]
  
  bioregion_data <- data[["bioreg_sf"]]
  colnames(eps_bioregion) <- bioregion_data$bioregion
  
  plot <- mcmc_areas(eps_bioregion, prob_outer = 0.9) +
    ggtitle(species)
  
  return(plot)
}

bio_plot_data <- list()

for(i in 1:length(species_names)) {
 bio_plot_data[[i]] <- bioregion_contribution(model = model_fits[[top]], 
                         data = multispecies_data, 
                         species = species_names[i], species_index = i)
} 

cowplot::plot_grid(plotlist = bio_plot_data, ncol = 1)
Influence of the bioregion random effect on abundance (log-scale).

Figure 12: Influence of the bioregion random effect on abundance (log-scale).

Bioregion usually had a substantial impact on abundance estimates across the state. The variance (standard deviation) associated with the bioregion random effect was large for Sambar (\(\sigma\) = 3.96, 90% CI: 2.8 - 5.38), but less for Fallow (\(\sigma\) = 1.67, 90% CI: 0.85 - 2.73). Red Deer also exhibited large variance in the bioregion random effect (\(\sigma\) = 3.7, 90% CI: 2.48 - 5.23), as well Hog (\(\sigma\) = 4.4, 90% CI: 3.11 - 5.93). The larger variances for Hog, Red and Sambar reflect their more refined range boundaries, wheras Fallow show less clear range boundaries and thus show less variance across the state.

Show code
# Bioregion map 
bioregion_eff <- model_fits[[top]]$summary("eps_bioregion", posterior::default_summary_measures())

bioreg <- VicmapR::vicmap_query("open-data-platform:vbioreg100") %>%
    VicmapR::collect() %>%
  group_by(bioregion, bioregion_code) %>%
  summarise(geometry = st_union(geometry)) %>%
  ungroup()

bioregion_nhat <- model_fits[[top]]$summary("Nhat_bioreg", prob_occ = ~ mean(. > 0), trimmed_mean = ~mean(., trim = 0.025), posterior::default_summary_measures())

bioregion_density <- list()

for(i in 1:length(species_names)) {
  which_sp <- which(stringr::str_detect(string = bioregion_nhat$variable,
                                         pattern = paste0("Nhat_bioreg\\[", i)))
  
  bioregion_nhat_sp <- bioregion_nhat[which_sp,]
  
  bioregion_nhat_sp$bioregion <- multispecies_data$bioreg_sf$bioregion 
  bioregion_nhat_sp$Area <- as.numeric(table(multispecies_data$pred_bioreg))
  
  bioregion_nhat_sp <- bioregion_nhat_sp %>% 
    mutate(`Average Density` = trimmed_mean/Area, 
           `Density (5%)` = q5/Area, 
           `Density (95%)` = q95/Area, 
           `CV` = sd/trimmed_mean, 
           label = paste0(bioregion, " [D = "
                          ,round(`Average Density`,2),
                          " (",round(`Density (5%)`,2)," - ",round(`Density (95%)`, 2),")]"))
  
  bioregion_density[[i]] <- left_join(bioreg, bioregion_nhat_sp, by = "bioregion")
  
}

bioregion_effects <- list()

for(i in 1:length(species_names)) {
  which_sp <- which(stringr::str_detect(string = bioregion_eff$variable,
                                         pattern = paste0("eps_bioregion\\[", i)))
  
  bioregion_eff_sp <- bioregion_eff[which_sp,]
  
  intercept <-  model_fits[[top]]$summary(paste0("beta_psi[",i,",1]"), mean = ~ mean(.))$mean
  
  bioregion_eff_sp$bioregion <- multispecies_data$bioreg_sf$bioregion 
  bioregion_eff_sp$Area <- as.numeric(table(multispecies_data$pred_bioreg))
  
  bioregion_eff_sp <- bioregion_eff_sp %>% 
    mutate(`Average Effect` = mean)
  
  bioregion_effects[[i]] <- left_join(bioreg, bioregion_eff_sp, by = "bioregion")
  
}

bioregion_ggplot <- function(effect_df, title) {
  ggplot(data = effect_df) +
    geom_sf(data = state, linewidth = 0.5, fill = "grey90") +
    geom_sf(aes(fill = `Average Effect`)) +
    geom_label_repel(data = effect_df, aes(label = bioregion, fill = `Average Effect`, geometry = geometry), 
                  size = 2.5, stat = "sf_coordinates", colour = "white") +
    scale_fill_gradient2(low = "#B9C600",
                         mid = delwp_cols[["Teal"]], 
                         high = delwp_cols[["Navy"]], guide = guide_coloursteps(title = "Average Effect\n (log-scale)")) +
    labs(x = "Longitude", y = "Latitude") +
    ggspatial::annotation_scale() +
    delwp_theme()
}

# cl_br <- st_intersection(crown_land, bioreg %>% filter(bioregion == "Bridgewater") %>% st_transform(3111))

cpb <- cowplot::plot_grid(bioregion_ggplot(bioregion_effects[[1]], "Sambar"),
bioregion_ggplot(bioregion_effects[[2]], "Fallow"),
bioregion_ggplot(bioregion_effects[[3]], "Red"),
bioregion_ggplot(bioregion_effects[[4]], "Hog"), labels = c("A) Sambar Deer", "B) Fallow Deer", "C) Red Deer", "D) Hog Deer"), ncol = 2)

ggsave(plot = cpb, filename = "outputs/plots/bioregion_effects_plot.png", width = 4800, height = 3200, units = "px", bg = "white")

cpb
Effect of Victorian Bioregions on abundance of deer

Figure 13: Effect of Victorian Bioregions on abundance of deer

Fixed-parameter Effects

Show code
# beta_summaries 
beta_summaries <- model_fits[[top]]$summary("beta_psi")

beta_table_summary <- beta_summaries %>% 
  mutate(species = factor(rep(species_names, times = 6),levels = species_names), 
         covariate = rep(c("(Intercept)",
              "Bare Soil (%)", 
              "Nitrogen (%)", 
              "Distance to Pastural Land (m)", 
              "Precipitation Seasonality", 
              'Forest Edge per km^2^ (m)'), each = 4)) %>%
  arrange(species) %>%
  select(species, covariate, mean, sd, q5, q95)

format_cov <- function(sp, b, data = beta_summaries, variable = "beta_psi") {
  fd <- data %>%
    filter(variable == !!paste0(variable, "[", sp, ",",b+1,"]"))
  
  paste0("$\\beta$ = ",
         round(fd$mean, 2), 
        " [90% CI: ", 
        round(fd$q5, 2), 
        " – ",
        round(fd$q95, 2),
        "]")
}

ab_joined_list <- list()
for(j in 1:length(species_names)) {
beta_draws <- model_fits[[top]]$draws("beta_psi", format = "df")

# which_sp <- which(stringr::str_detect(string = colnames(beta_draws),
#                                     pattern = paste0("beta_psi\\[", i)))
# 
# beta_draws <- beta_draws[,which_sp]

ab_marginal_effects <- list()
ab_plot <- list()

ab_to_use <- ab_formula_3

phi_vars <- unlist(stringr::str_remove_all(string = labels(terms(ab_to_use)), pattern =  "log|scale|sqrt|\\)|\\("))
phi_logs <- unlist(stringr::str_detect(string = labels(terms(ab_to_use)), pattern =  "log\\("))
phi_sqrts <- c(0.5,0.5,1,1,0.5)

phi_labs <- c("Bare Soil (%)", 
              "Nitrogen (%)", 
              "Distance to Pastural Land (m)", 
              "Precipitation Seasonality", 
              'Forest Edge per km^2^ (m)')

fac <- c(FALSE, FALSE, FALSE, FALSE, FALSE)

for(i in 1:length(phi_vars)) {
ab_marginal_effects[[i]] <- marginal_effects_cmd(beta_draws, 
                                              param = "beta_psi", species_index = j,
                                              param_number = i+1, log = phi_logs[i],
                                     model_data = multispecies_data[["raw_data"]] %>% mutate(Nitrogen = Nitrogen*100), 
                                     abundance = TRUE,
                                     pwr = phi_sqrts[i],
                                     model_column = phi_vars[i], 
                                     transition = FALSE) %>%
  mutate(species = species_names[j])

}
ab_joined_list[[j]] <- bind_rows(ab_marginal_effects)
}

all_me_data <- bind_rows(ab_joined_list) %>%
  mutate(param = case_when(param == "BareSoil" ~ "Bare Soil (%)", 
                           param == "Nitrogen" ~ "Nitrogen (%)", 
                           param == "PastureDistance" ~ "Distance to Pastural Land (m)", 
                           param == "BIO15" ~ "Precipitation Seasonality", 
                           param == "ForestEdge" ~ 'Forest Edge per km^2^ (m)'))
  

cov_plot <- marginal_effects_plot_cmd_all(all_me_data %>% mutate(species = factor(species, levels = species_names), 
                                                     param = factor(param, levels = phi_labs)), 
                          col = "DarkGreen", 
                          factor = FALSE,
                          ylab = "Contribution to Abundance (log-scale)") +
      ggplot2::facet_grid(species~param, scales = "free") +
  delwp_theme() +
  theme(strip.text = ggtext::element_markdown()) +
  xlab("Covariate Value")

ggsave(plot = cov_plot, filename = "outputs/plots/abundance_covariate_effects.png", width = 3000, height = 2000, units = "px")

cov_plot
Marginal response curves for the various fixed-effect parameters used in the model.

Figure 14: Marginal response curves for the various fixed-effect parameters used in the model.

Of the five fixed-effect covariates used in the top model we found that all had tangible effects for one or more species. The direction and magnitude of these effects often varied between species. Bare soil had a slight positive effect on Sambar abundance (\(\beta\) = 0.37 [90% CI: 0.02 – 0.72]), negligible effect on Fallow (\(\beta\) = -0.1 [90% CI: -0.63 – 0.4]) and Red Deer (\(\beta\) = -0.73 [90% CI: -2.01 – 0.54]), and a slight negative effect on Hog Deer abundance (\(\beta\) = -0.72 [90% CI: -1.57 – 0.12]). Nitrogen generally had weaker, more variable effects. The effect of nitrogen was weakly positive for Sambar (\(\beta\) = 0.31 [90% CI: -0.08 – 0.68]), and Hog Deer (\(\beta\) = 0.53 [90% CI: -0.54 – 1.6]); but weakly negative for Fallow (\(\beta\) = 0.22 [90% CI: -0.36 – 0.72]), and Red Deer (\(\beta\) = -0.45 [90% CI: -1.48 – 0.63]). The effect of distance to pasture was consistently negative across all four species, which can also be interpreted as a prefence of deer to be ‘closer’ to pastural areas. However the effect of this covariate varied in magnitude for Sambar (\(\beta\) = -0.53 [90% CI: -0.83 – -0.22]), Fallow (\(\beta\) = -0.89 [90% CI: -1.35 – -0.43]), Red (\(\beta\) = -1.61 [90% CI: -2.71 – -0.48]), and Hog Deer (\(\beta\) = -0.3 [90% CI: -0.8 – 0.22]). Sambar Deer weakly favored areas with less seasonal variation in precipitation (\(\beta\) = -0.23 [90% CI: -0.57 – 0.12]), with Fallow (\(\beta\) = 0.39 [90% CI: 0.02 – 0.8]), Red (\(\beta\) = 1.05 [90% CI: 0.05 – 2.02]), and Hog Deer (\(\beta\) = 0.77 [90% CI: 0.21 – 1.31]) all showing positive relationships between precipitation seasonality and abundance. Lastly, we saw that the amount of forest edge in the landscape positively impacted abundance for Sambar (\(\beta\) = 0.17 [90% CI: -0.06 – 0.4]), and Fallow Deer (\(\beta\) = 0.61 [90% CI: 0.29 – 0.93]). This relationship was weakly negative for Red Deer (\(\beta\) = -0.7 [90% CI: -1.39 – 0.01]), and weakly positive for Hog Deer (\(\beta\) = 0.32 [90% CI: -0.07 – 0.69]).

Show code
which_sig <- which(beta_table_summary$q5 * beta_table_summary$q95 > 0)

write.csv(beta_table_summary, "outputs/tables/beta_coefs.csv")
beta_table_summary %>%
  kbl("html", digits = 3, caption = "beta-coeffiient estimates for the drivers of species abundance, coefficient estimates not-overalpping zero are shown in bold") %>%
  kable_styling("striped") %>%
  collapse_rows(1) %>%
  row_spec(row = which_sig, bold = TRUE)
Table 6: beta-coeffiient estimates for the drivers of species abundance, coefficient estimates not-overalpping zero are shown in bold
species covariate mean sd q5 q95
Sambar Deer (Intercept) -3.615 0.956 -5.182 -2.166
Bare Soil (%) 0.368 0.215 0.023 0.723
Nitrogen (%) 0.313 0.230 -0.079 0.683
Distance to Pastural Land (m) -0.535 0.189 -0.832 -0.221
Precipitation Seasonality -0.230 0.208 -0.573 0.122
Forest Edge per km2 (m) 0.166 0.137 -0.057 0.396
Fallow Deer (Intercept) -1.680 0.513 -2.563 -0.930
Bare Soil (%) -0.103 0.319 -0.630 0.403
Nitrogen (%) 0.220 0.335 -0.359 0.723
Distance to Pastural Land (m) -0.892 0.278 -1.345 -0.432
Precipitation Seasonality 0.392 0.239 0.023 0.799
Forest Edge per km2 (m) 0.607 0.192 0.288 0.931
Red Deer (Intercept) -6.781 1.170 -8.861 -4.997
Bare Soil (%) -0.725 0.793 -2.013 0.543
Nitrogen (%) -0.447 0.635 -1.481 0.628
Distance to Pastural Land (m) -1.613 0.683 -2.705 -0.482
Precipitation Seasonality 1.049 0.599 0.052 2.021
Forest Edge per km2 (m) -0.699 0.426 -1.386 0.005
Hog Deer (Intercept) -8.297 1.307 -10.503 -6.305
Bare Soil (%) -0.717 0.512 -1.567 0.120
Nitrogen (%) 0.526 0.651 -0.543 1.601
Distance to Pastural Land (m) -0.298 0.312 -0.799 0.222
Precipitation Seasonality 0.772 0.334 0.209 1.312
Forest Edge per km2 (m) 0.318 0.231 -0.068 0.690

Additional Parameters of Interest

Below is a table with several model parameters we estimated:

Show code
key_params <- model_fits[[top]]$summary(c("av_gs", 
                                          "od_mu", 
                                          "odRN", 
                                          "activ", 
                                          "Nhat_sum"))

key_params %>% 
  kbl(format = "html", digits = 2, caption = "Several key outputs and parameters of the model not reported earlier") %>%
  kable_styling("striped")
Table 7: Several key outputs and parameters of the model not reported earlier
variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
av_gs[1] 1.01 1.00 0.00 0.00 1.00 1.01 1.00 1067.27 1375.61
av_gs[2] 1.02 1.02 0.02 0.01 1.00 1.05 1.02 564.16 775.02
av_gs[3] 1.04 1.03 0.04 0.02 1.01 1.10 1.00 1923.67 1447.03
av_gs[4] 1.03 1.01 0.04 0.01 1.00 1.10 1.00 591.52 583.61
od_mu[1] -1.93 -1.93 0.21 0.21 -2.28 -1.59 1.00 965.07 1255.70
od_mu[2] -3.62 -3.61 0.27 0.26 -4.07 -3.18 1.01 774.65 1157.38
od_mu[3] -2.86 -2.86 0.36 0.37 -3.43 -2.29 1.01 1207.00 1323.88
od_mu[4] -2.54 -2.53 0.32 0.33 -3.08 -2.00 1.00 913.74 916.78
odRN 0.79 0.78 0.14 0.14 0.59 1.04 1.01 815.79 1125.54
activ 0.52 0.52 0.02 0.02 0.48 0.56 1.00 2217.73 1339.29
Nhat_sum 192366.95 186810.00 34001.78 30296.93 146731.70 255490.00 1.00 1450.13 1403.75

Dispersion

Below we show the parameter estimates for the dispersion parameters estimated within our negative binomial models for camera counts and Royle-Nichols model.

Show code
dispersion_data <- model_fits[[top]]$summary(c("od", "odRN")) %>%
  mutate(species = c(species_names, "All"), 
         model = c(rep("Camera count", 4), "Royle-Nichols"), 
         parameter = c("$\\phi_1$", "$\\phi_2$", "$\\phi_3$", "$\\phi_4$", "$\\bar \\phi$")) %>%
  select(model, species, parameter, mean, sd, q5, q95)

format_od <- function(data) {
  paste0(dispersion_data$parameter, " = ", round(dispersion_data$mean, 3), " [90 % CI: ", round(dispersion_data$q5,3), " -- ", round(dispersion_data$q95,3), "]")
}

formatred_od <- format_od(dispersion_data)

dispersion_data %>% 
  kbl("html", digits = 3, caption = "Estimates of the dispersion parameter within negative binomial submodels") %>%
  kable_styling("striped") %>%
  collapse_rows(1)
Table 8: Estimates of the dispersion parameter within negative binomial submodels
model species parameter mean sd q5 q95
Camera count Sambar Deer \(\phi_1\) 0.114 0.014 0.093 0.137
Fallow Deer \(\phi_2\) 0.021 0.004 0.015 0.030
Red Deer \(\phi_3\) 0.047 0.016 0.025 0.076
Hog Deer \(\phi_4\) 0.064 0.018 0.038 0.095
Royle-Nichols All \(\bar \phi\) 0.790 0.140 0.591 1.041

Estimates of abundance were also dependent upon overdispersion (excess variation) in the camera counts and transect detections. Sambar (\(\phi_1\) = 0.114 [90 % CI: 0.093 – 0.137]), Fallow (\(\phi_2\) = 0.021 [90 % CI: 0.015 – 0.03]), Red (\(\phi_3\) = 0.047 [90 % CI: 0.025 – 0.076]), and Hog Deer (\(\phi_4\) = 0.064 [90 % CI: 0.038 – 0.095]) all had large amounts of overdispersion (where lower values (<1) equate to more variance). The Royle-Nichols model had less variance (than camera counts) associated with transect-level detections/non-detections (\(\bar \phi\) = 0.79 [90 % CI: 0.591 – 1.041]).

Relating abundance to vegetation

From our point density estimates of deer, we can model relationships between deer densities at sites and a range of structural vegetation measures collected during the surveys. For this part of the analysis we implement a multivariate regression.

Firstly the response variables we model are:

While our main predictor of interest is the density of deer at a site, we also control for predictors that are expected to have an impact on vegetation growth and composition. The predictors controlled for include features relating to climate, ecology, topography, fire history and forest structure. They are the following predictors:

Data Preparation

We can use the sum of estimated deer density for the four species at each site as our key predictor variable. Given that we already extracted a range of bioclimatic and environmental covariates for the original model, we can combine those with the data available here.

Note that this model is run on a subset of sites (sites surveyed as part of the statewide deer survey), as the hog deer sites did not have complete vegetation metrics recorded.

Show code
density <- model_fits[[top]]$draws("N_site", format = "matrix")
site_estimates <- list()
site_density <- list()

for(i in 1:n_site) {
  site_estimates[[i]] <- character()
  for(j in 1:length(deer_species_all)) {
    site_estimates[[i]][j] <- paste0("N_site[", j, ",", i, "]")
  }
  site_density[[i]] <- sum(colMeans(density)[site_estimates[[i]]])
}

site_density_ul <- unlist(site_density)

site_density_df <- cams_curated %>%
  mutate(AllDeerDensity = site_density_ul, 
         Bioregion = factor(multispecies_data$site_bioreg), 
         EVC = factor(multispecies_data$site_evc)) %>% 
  filter(ProjectShortName == "StatewideDeer") %>% 
  left_join(multispecies_data[["raw_data"]] %>% select(-EVC, - BIOREGION)) %>%
  mutate(TSLF_bin = as.factor(round(TSLF_bin)))

 site_vars <- dplyr::tbl(con, dbplyr::in_schema("deervic", "curated_site_data")) %>%
    dplyr::filter(SiteID %in% !!c(cams_curated$SiteID, "47191")) %>%
    dplyr::collect() %>%
    dplyr::mutate(HerbaceousUnderstoryCover = NNWHUCover + ENWHUCover,
                  Exotics = EWUCover + ENWHUCover,
                  ExoticPresence = as.integer(as.logical(Exotics)),
                  SeedlingsOrSaplings = as.integer(as.logical(Seedlings+Saplings)),
                  SeedlingsPresence = as.integer(as.logical(Seedlings)),
                  SaplingsPresence = as.integer(as.logical(Saplings)),
                  SiteID = dplyr::case_when(SiteID == "47191" & CameraID == "HO04101053" ~ "47191A",
                                            SiteID == "47191" & CameraID != "HO04101053" ~ "47191B",
                                            TRUE ~ SiteID)) %>% # native + exotic herbaceous cover
    dplyr::arrange(SiteID) %>%
    as.data.frame()
 
 site_density_df_joined <- site_density_df %>% 
   left_join(site_vars) %>%
   mutate(BGroundCover = ifelse(BGroundCover == 0, 0.125, BGroundCover)/100, 
             NWUCover = ifelse(NWUCover == 0, 0.125, NWUCover)/100, 
             NNWHUCover = ifelse(NNWHUCover == 0, 0.125, NNWHUCover)/100)

Model Execution

We use the brms package to run our model in a Bayesian framework (Bürkner 2017, 2018, 2021). Our model is run in parallel on four cores over four chains for 1,000 sampling iterations per chain (1,000 warm up iterations).

Show code
bform1 <- bf(mvbind(BGroundCover, 
                    NWUCover, 
                    NNWHUCover, 
                    Seedlings,
                    Saplings,
                    ExoticPresence) ~ AllDeerDensity
                                    + scale(BIO01)
                                    + scale(BIO04)
                                    + scale(BIO06)
                                    + scale(BIO12)
                                    + scale(BIO15)
                                    + scale(Nitrogen)
                                    + TSLF_bin
                                    + scale(TWIND)
                                    + scale(sqrt(ForestEdge))
                                    + scale(sqrt(CanopyCov))
                                    + scale(sqrt(TopHeight))
                                    + (1|EVC)) 


vegfit <- brm(bform1, 
            data = site_density_df_joined, 
            family = list(Beta(link = "logit"), 
                          Beta(link = "logit"),
                          Beta(link = "logit"),
                          negbinomial(),
                          negbinomial(),
                          bernoulli(link = "logit")),
            inits = "0",
            chains = 4, 
            cores = 4, 
            iter = 2000, 
            warmup = 1000,
            backend = "cmdstanr")

saveRDS(vegfit, "outputs/models/vegfit.rds")

Model Results

Model Summary

Below we show a complete summary of the model, including all coefficients estimated for all the parameters and response variables. Standard deviations for the EVC random effect and the estimates of the \(\phi\) parameter in the beta-distribution models, and the shape parameter in the negative binomial models are also shown. Note that all Rhats > 1.05 and all ESS > 200, suggesting good chain mixing and convergence.

Show code
vegfit <- readRDS("outputs/models/vegfit.rds")
veg_coefs <- fixef(vegfit, probs = c(0.05, 0.95), pars = c("BGroundCover_AllDeerDensity", 
                                                             "NWUCover_AllDeerDensity", 
                                                             "NNWHUCover_AllDeerDensity", 
                                                             "Seedlings_AllDeerDensity", 
                                                             "Saplings_AllDeerDensity", 
                                                             "ExoticPresence_AllDeerDensity")) %>%
  as.data.frame() %>%
  mutate(Effect = paste0(c("Bare Ground Cover",
                           "Native Woody Understorey Cover",
                           "Native Non-woody Understorey Cover", 
                           "Presence of seedlings",
                           "Presence of saplings",
                           "Presence of exotic flora"), " ~ Deer Density"))

print_conditional <- function(x, i, d = 3, OR = FALSE) {
  row <- x[i,]
  name <- "\\beta"
  if(OR) {
    row <- row %>%
      mutate(across(where(is.numeric), exp))
    name <- "OR"
  }
  
paste0("$",name,"$ = ", round(row[["Estimate"]], d), " [90% CI: ", round(row[["Q5"]], d), ", ", round(row[["Q95"]], d), "]")
}

veg_coefs_save <- veg_coefs %>%
  transmute(effect = Effect, mean = Estimate, sd = Est.Error, `5 %` = Q5, `95 %` = Q95) %>%
  mutate(across(where(is.numeric), round, 3))

write.csv(veg_coefs_save, "outputs/tables/veg_coefs.csv", row.names = F)


summary(vegfit) 
 Family: MV(beta, beta, beta, negbinomial, negbinomial, bernoulli) 
  Links: mu = logit; phi = identity
         mu = logit; phi = identity
         mu = logit; phi = identity
         mu = log; shape = identity
         mu = log; shape = identity
         mu = logit 
Formula: BGroundCover ~ AllDeerDensity + scale(BIO01) + scale(BIO04) + scale(BIO06) + scale(BIO12) + scale(BIO15) + scale(Nitrogen) + TSLF_bin + scale(TWIND) + scale(sqrt(ForestEdge)) + scale(sqrt(CanopyCov)) + scale(sqrt(TopHeight)) + (1 | EVC) 
         NWUCover ~ AllDeerDensity + scale(BIO01) + scale(BIO04) + scale(BIO06) + scale(BIO12) + scale(BIO15) + scale(Nitrogen) + TSLF_bin + scale(TWIND) + scale(sqrt(ForestEdge)) + scale(sqrt(CanopyCov)) + scale(sqrt(TopHeight)) + (1 | EVC) 
         NNWHUCover ~ AllDeerDensity + scale(BIO01) + scale(BIO04) + scale(BIO06) + scale(BIO12) + scale(BIO15) + scale(Nitrogen) + TSLF_bin + scale(TWIND) + scale(sqrt(ForestEdge)) + scale(sqrt(CanopyCov)) + scale(sqrt(TopHeight)) + (1 | EVC) 
         Seedlings ~ AllDeerDensity + scale(BIO01) + scale(BIO04) + scale(BIO06) + scale(BIO12) + scale(BIO15) + scale(Nitrogen) + TSLF_bin + scale(TWIND) + scale(sqrt(ForestEdge)) + scale(sqrt(CanopyCov)) + scale(sqrt(TopHeight)) + (1 | EVC) 
         Saplings ~ AllDeerDensity + scale(BIO01) + scale(BIO04) + scale(BIO06) + scale(BIO12) + scale(BIO15) + scale(Nitrogen) + TSLF_bin + scale(TWIND) + scale(sqrt(ForestEdge)) + scale(sqrt(CanopyCov)) + scale(sqrt(TopHeight)) + (1 | EVC) 
         ExoticPresence ~ AllDeerDensity + scale(BIO01) + scale(BIO04) + scale(BIO06) + scale(BIO12) + scale(BIO15) + scale(Nitrogen) + TSLF_bin + scale(TWIND) + scale(sqrt(ForestEdge)) + scale(sqrt(CanopyCov)) + scale(sqrt(TopHeight)) + (1 | EVC) 
   Data: site_density_df_joined (Number of observations: 252) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Group-Level Effects: 
~EVC (Number of levels: 19) 
                             Estimate Est.Error l-95% CI u-95% CI
sd(BGroundCover_Intercept)       0.28      0.11     0.07     0.54
sd(NWUCover_Intercept)           0.14      0.10     0.01     0.37
sd(NNWHUCover_Intercept)         0.22      0.13     0.01     0.51
sd(Seedlings_Intercept)          1.18      0.40     0.51     2.09
sd(Saplings_Intercept)           0.84      0.33     0.17     1.56
sd(ExoticPresence_Intercept)     0.55      0.38     0.03     1.45
                             Rhat Bulk_ESS Tail_ESS
sd(BGroundCover_Intercept)   1.01     1208      924
sd(NWUCover_Intercept)       1.00     1542     2250
sd(NNWHUCover_Intercept)     1.01     1143     1630
sd(Seedlings_Intercept)      1.00     1138     1655
sd(Saplings_Intercept)       1.01      727      643
sd(ExoticPresence_Intercept) 1.00      724     1352

Population-Level Effects: 
                                   Estimate Est.Error l-95% CI
BGroundCover_Intercept                -2.13      0.22    -2.58
NWUCover_Intercept                    -0.53      0.18    -0.87
NNWHUCover_Intercept                  -0.47      0.22    -0.91
Seedlings_Intercept                    3.26      0.58     2.11
Saplings_Intercept                     3.50      0.49     2.50
ExoticPresence_Intercept              -1.72      0.55    -2.83
BGroundCover_AllDeerDensity           -0.04      0.03    -0.10
BGroundCover_scaleBIO01                0.37      0.42    -0.46
BGroundCover_scaleBIO04                0.08      0.26    -0.43
BGroundCover_scaleBIO06               -0.22      0.48    -1.15
BGroundCover_scaleBIO12               -0.27      0.15    -0.55
BGroundCover_scaleBIO15               -0.20      0.09    -0.36
BGroundCover_scaleNitrogen            -0.01      0.16    -0.32
BGroundCover_TSLF_bin3                 0.16      0.22    -0.26
BGroundCover_TSLF_bin4                 0.11      0.24    -0.36
BGroundCover_TSLF_bin5                -0.19      0.24    -0.66
BGroundCover_scaleTWIND                0.06      0.11    -0.15
BGroundCover_scalesqrtForestEdge      -0.04      0.07    -0.18
BGroundCover_scalesqrtCanopyCov       -0.02      0.07    -0.15
BGroundCover_scalesqrtTopHeight       -0.26      0.08    -0.41
NWUCover_AllDeerDensity               -0.05      0.03    -0.10
NWUCover_scaleBIO01                    0.08      0.41    -0.74
NWUCover_scaleBIO04                   -0.04      0.26    -0.54
NWUCover_scaleBIO06                   -0.08      0.48    -1.00
NWUCover_scaleBIO12                    0.18      0.14    -0.10
NWUCover_scaleBIO15                    0.19      0.08     0.03
NWUCover_scaleNitrogen                 0.02      0.16    -0.30
NWUCover_TSLF_bin3                    -0.35      0.19    -0.73
NWUCover_TSLF_bin4                    -0.58      0.22    -1.02
NWUCover_TSLF_bin5                    -0.91      0.24    -1.37
NWUCover_scaleTWIND                   -0.01      0.12    -0.24
NWUCover_scalesqrtForestEdge           0.05      0.08    -0.10
NWUCover_scalesqrtCanopyCov           -0.05      0.07    -0.19
NWUCover_scalesqrtTopHeight           -0.14      0.08    -0.29
NNWHUCover_AllDeerDensity              0.07      0.03     0.01
NNWHUCover_scaleBIO01                  0.25      0.48    -0.71
NNWHUCover_scaleBIO04                 -0.36      0.29    -0.92
NNWHUCover_scaleBIO06                 -0.58      0.55    -1.64
NNWHUCover_scaleBIO12                  0.42      0.16     0.09
NNWHUCover_scaleBIO15                  0.11      0.09    -0.07
NNWHUCover_scaleNitrogen              -0.31      0.17    -0.65
NNWHUCover_TSLF_bin3                  -0.32      0.23    -0.76
NNWHUCover_TSLF_bin4                  -0.15      0.25    -0.64
NNWHUCover_TSLF_bin5                  -0.15      0.27    -0.66
NNWHUCover_scaleTWIND                 -0.10      0.13    -0.36
NNWHUCover_scalesqrtForestEdge         0.01      0.08    -0.16
NNWHUCover_scalesqrtCanopyCov          0.04      0.08    -0.11
NNWHUCover_scalesqrtTopHeight          0.03      0.08    -0.13
Seedlings_AllDeerDensity               0.06      0.07    -0.08
Seedlings_scaleBIO01                   0.15      0.92    -1.67
Seedlings_scaleBIO04                  -0.12      0.58    -1.27
Seedlings_scaleBIO06                  -0.58      1.01    -2.55
Seedlings_scaleBIO12                   0.67      0.43    -0.11
Seedlings_scaleBIO15                  -0.22      0.28    -0.78
Seedlings_scaleNitrogen               -0.92      0.40    -1.71
Seedlings_TSLF_bin3                   -1.19      0.53    -2.25
Seedlings_TSLF_bin4                   -1.22      0.55    -2.34
Seedlings_TSLF_bin5                   -1.64      0.59    -2.81
Seedlings_scaleTWIND                  -0.07      0.24    -0.52
Seedlings_scalesqrtForestEdge         -0.36      0.18    -0.70
Seedlings_scalesqrtCanopyCov          -0.09      0.18    -0.44
Seedlings_scalesqrtTopHeight          -0.90      0.21    -1.31
Saplings_AllDeerDensity               -0.01      0.06    -0.13
Saplings_scaleBIO01                    0.28      0.87    -1.35
Saplings_scaleBIO04                    0.08      0.53    -0.96
Saplings_scaleBIO06                    0.16      0.95    -1.74
Saplings_scaleBIO12                    0.29      0.32    -0.34
Saplings_scaleBIO15                   -0.42      0.19    -0.78
Saplings_scaleNitrogen                 0.14      0.34    -0.51
Saplings_TSLF_bin3                    -0.56      0.48    -1.47
Saplings_TSLF_bin4                    -1.76      0.49    -2.72
Saplings_TSLF_bin5                    -1.30      0.54    -2.32
Saplings_scaleTWIND                    0.09      0.24    -0.36
Saplings_scalesqrtForestEdge          -0.31      0.16    -0.62
Saplings_scalesqrtCanopyCov            0.02      0.14    -0.26
Saplings_scalesqrtTopHeight           -0.56      0.17    -0.90
ExoticPresence_AllDeerDensity          0.16      0.07     0.02
ExoticPresence_scaleBIO01              0.34      1.11    -1.81
ExoticPresence_scaleBIO04             -0.11      0.71    -1.52
ExoticPresence_scaleBIO06             -0.65      1.30    -3.24
ExoticPresence_scaleBIO12              0.37      0.39    -0.38
ExoticPresence_scaleBIO15              0.34      0.21    -0.07
ExoticPresence_scaleNitrogen          -1.29      0.46    -2.23
ExoticPresence_TSLF_bin3               0.55      0.58    -0.62
ExoticPresence_TSLF_bin4               0.64      0.64    -0.62
ExoticPresence_TSLF_bin5               0.95      0.65    -0.37
ExoticPresence_scaleTWIND             -0.45      0.30    -1.07
ExoticPresence_scalesqrtForestEdge     0.14      0.19    -0.22
ExoticPresence_scalesqrtCanopyCov      0.27      0.19    -0.10
ExoticPresence_scalesqrtTopHeight      0.12      0.21    -0.29
                                   u-95% CI Rhat Bulk_ESS Tail_ESS
BGroundCover_Intercept                -1.71 1.00     2209     3079
NWUCover_Intercept                    -0.18 1.00     3213     2770
NNWHUCover_Intercept                  -0.04 1.00     2455     2646
Seedlings_Intercept                    4.40 1.00     2090     2639
Saplings_Intercept                     4.43 1.00     2097     2186
ExoticPresence_Intercept              -0.66 1.00     2797     2798
BGroundCover_AllDeerDensity            0.01 1.00     5036     3211
BGroundCover_scaleBIO01                1.17 1.00     2277     2814
BGroundCover_scaleBIO04                0.60 1.00     2248     2990
BGroundCover_scaleBIO06                0.71 1.00     2310     2772
BGroundCover_scaleBIO12                0.02 1.00     4669     3069
BGroundCover_scaleBIO15               -0.04 1.00     3370     3575
BGroundCover_scaleNitrogen             0.31 1.00     4330     2948
BGroundCover_TSLF_bin3                 0.59 1.00     2406     3050
BGroundCover_TSLF_bin4                 0.58 1.00     2493     3072
BGroundCover_TSLF_bin5                 0.29 1.00     2871     2940
BGroundCover_scaleTWIND                0.28 1.00     5594     3415
BGroundCover_scalesqrtForestEdge       0.11 1.00     5232     3414
BGroundCover_scalesqrtCanopyCov        0.12 1.01     5684     3341
BGroundCover_scalesqrtTopHeight       -0.10 1.00     6654     3086
NWUCover_AllDeerDensity                0.00 1.00     5754     2689
NWUCover_scaleBIO01                    0.88 1.00     1832     2781
NWUCover_scaleBIO04                    0.47 1.00     1885     2862
NWUCover_scaleBIO06                    0.87 1.00     1848     2817
NWUCover_scaleBIO12                    0.47 1.00     3838     3044
NWUCover_scaleBIO15                    0.36 1.00     3812     3106
NWUCover_scaleNitrogen                 0.34 1.00     4381     3131
NWUCover_TSLF_bin3                     0.02 1.00     3273     3393
NWUCover_TSLF_bin4                    -0.13 1.00     3033     3305
NWUCover_TSLF_bin5                    -0.45 1.00     3032     3043
NWUCover_scaleTWIND                    0.22 1.00     5076     2759
NWUCover_scalesqrtForestEdge           0.19 1.00     6283     3518
NWUCover_scalesqrtCanopyCov            0.09 1.00     5756     3482
NWUCover_scalesqrtTopHeight            0.01 1.00     5135     3027
NNWHUCover_AllDeerDensity              0.12 1.00     5249     3097
NNWHUCover_scaleBIO01                  1.17 1.00     1660     2327
NNWHUCover_scaleBIO04                  0.22 1.00     1742     2245
NNWHUCover_scaleBIO06                  0.48 1.00     1684     2216
NNWHUCover_scaleBIO12                  0.73 1.00     4083     3334
NNWHUCover_scaleBIO15                  0.27 1.00     3854     2788
NNWHUCover_scaleNitrogen               0.02 1.00     4044     3426
NNWHUCover_TSLF_bin3                   0.12 1.00     2698     2614
NNWHUCover_TSLF_bin4                   0.36 1.00     2491     2805
NNWHUCover_TSLF_bin5                   0.39 1.00     2644     2777
NNWHUCover_scaleTWIND                  0.16 1.00     4115     2774
NNWHUCover_scalesqrtForestEdge         0.17 1.00     5317     2833
NNWHUCover_scalesqrtCanopyCov          0.19 1.00     5677     3427
NNWHUCover_scalesqrtTopHeight          0.19 1.00     4446     2908
Seedlings_AllDeerDensity               0.20 1.00     4816     3274
Seedlings_scaleBIO01                   1.92 1.00     1956     2536
Seedlings_scaleBIO04                   1.04 1.00     2077     2446
Seedlings_scaleBIO06                   1.43 1.00     1989     2647
Seedlings_scaleBIO12                   1.50 1.00     2099     2597
Seedlings_scaleBIO15                   0.34 1.00     2835     2727
Seedlings_scaleNitrogen               -0.13 1.00     4481     2949
Seedlings_TSLF_bin3                   -0.16 1.00     2719     3097
Seedlings_TSLF_bin4                   -0.16 1.00     2345     2680
Seedlings_TSLF_bin5                   -0.48 1.00     2085     2893
Seedlings_scaleTWIND                   0.41 1.00     3846     3337
Seedlings_scalesqrtForestEdge         -0.01 1.00     4923     2917
Seedlings_scalesqrtCanopyCov           0.27 1.00     4579     2934
Seedlings_scalesqrtTopHeight          -0.51 1.00     5527     3278
Saplings_AllDeerDensity                0.10 1.00     5183     3482
Saplings_scaleBIO01                    2.00 1.00     1605     1905
Saplings_scaleBIO04                    1.09 1.00     1666     1889
Saplings_scaleBIO06                    1.95 1.00     1644     1915
Saplings_scaleBIO12                    0.90 1.00     2701     2994
Saplings_scaleBIO15                   -0.05 1.00     3362     2985
Saplings_scaleNitrogen                 0.80 1.00     3934     3542
Saplings_TSLF_bin3                     0.38 1.00     2534     3119
Saplings_TSLF_bin4                    -0.75 1.00     2467     2689
Saplings_TSLF_bin5                    -0.22 1.00     2400     3064
Saplings_scaleTWIND                    0.57 1.00     3318     3347
Saplings_scalesqrtForestEdge          -0.00 1.00     5403     3336
Saplings_scalesqrtCanopyCov            0.30 1.00     4816     2719
Saplings_scalesqrtTopHeight           -0.24 1.00     3151     2805
ExoticPresence_AllDeerDensity          0.31 1.00     4510     3108
ExoticPresence_scaleBIO01              2.54 1.00     2173     2623
ExoticPresence_scaleBIO04              1.27 1.00     2136     2708
ExoticPresence_scaleBIO06              1.84 1.00     2106     2689
ExoticPresence_scaleBIO12              1.16 1.00     4134     2957
ExoticPresence_scaleBIO15              0.75 1.00     4024     3078
ExoticPresence_scaleNitrogen          -0.42 1.00     3880     3239
ExoticPresence_TSLF_bin3               1.70 1.00     2252     3030
ExoticPresence_TSLF_bin4               1.86 1.00     2080     2968
ExoticPresence_TSLF_bin5               2.19 1.00     2240     2838
ExoticPresence_scaleTWIND              0.15 1.00     5481     3137
ExoticPresence_scalesqrtForestEdge     0.51 1.00     6163     3365
ExoticPresence_scalesqrtCanopyCov      0.65 1.00     5730     2985
ExoticPresence_scalesqrtTopHeight      0.52 1.00     5427     2947

Family Specific Parameters: 
                 Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
phi_BGroundCover     6.79      0.70     5.49     8.24 1.00     3784
phi_NWUCover         4.07      0.36     3.43     4.82 1.00     6565
phi_NNWHUCover       3.05      0.24     2.59     3.55 1.00     5868
shape_Seedlings      0.29      0.03     0.23     0.36 1.00     3508
shape_Saplings       0.38      0.04     0.30     0.46 1.00     3338
                 Tail_ESS
phi_BGroundCover     2773
phi_NWUCover         3006
phi_NNWHUCover       2838
shape_Seedlings      2573
shape_Saplings       2628

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Deer Density Conditional Effects

Given our main interest in the model is deer density and how it impacts the suite of response variables, we generate marginal/conditional effects plots for the relationship between deer and the vegetation metric. Based on these plots we can see a “significant” (confidence intervals not-overlapping zero) effect of deer density on: bare ground (-ve), native woody understorey (-ve), native non-woody/herbaceous understorey (+ve), and presence of exotic species (+ve):

Show code
plot_list <- conditional_effects(vegfit, effects = "AllDeerDensity", prob = 0.9)
plot_list_50 <- conditional_effects(vegfit, effects = "AllDeerDensity", prob = 0.5)

plot_deer_relationship <- function(plot_data, plot_data_50, yaxis) {
  ggplot(plot_data) +
    geom_ribbon(aes(x = AllDeerDensity, ymin = `lower__`, ymax = `upper__`), 
                fill = delwp_cols[["Navy"]], alpha = 0.3) +
    geom_ribbon(aes(x = AllDeerDensity, ymin = `lower__`, ymax = `upper__`), 
                fill = delwp_cols[["Navy"]], alpha = 0.3, data = plot_data_50) +
    geom_smooth(aes(x = AllDeerDensity, y = `estimate__`), 
                colour = delwp_cols[["Navy"]], method = "loess", formula = y~x, linewidth = 2) +
    ylab(paste(yaxis)) +
    xlab("Deer Density (per km^2^)") +
    delwp_theme() +
    theme(axis.title.x = ggtext::element_markdown())
    
}

vegplot <- cowplot::plot_grid(plot_deer_relationship(plot_list[[1]], plot_list_50[[1]], "Bare Ground Cover"), 
                   plot_deer_relationship(plot_list[[2]], plot_list_50[[2]], "Native Woody Understorey Cover"),
                   plot_deer_relationship(plot_list[[3]], plot_list_50[[3]], "Native Herbaceous Understorey Cover"),
                   plot_deer_relationship(plot_list[[4]], plot_list_50[[4]], "Seedling Counts") + scale_y_sqrt(),
                   plot_deer_relationship(plot_list[[5]], plot_list_50[[5]], "Sapling Counts"),
                   plot_deer_relationship(plot_list[[6]], plot_list_50[[6]], "Presence of Exotic Flora"), 
                   ncol = 3, labels = "AUTO")

ggsave(plot = vegplot, filename = "outputs/plots/vegetation_effects_plot.png", width = 3400, height = 2200, units = "px")

vegplot
Impact of deer density on a range of vegetation measure, including (A) bare ground cover, (B) native woody understorey, (C) native herbaceous understorey cover, (D) seedlings, (E) saplings, and (F) th presence of exotic species. The line shows the mean estimate, with the shading representing 50% and 90% confidence intervals. Inner and outer shading represents 50% and 90% confidence levels respectively.

Figure 15: Impact of deer density on a range of vegetation measure, including (A) bare ground cover, (B) native woody understorey, (C) native herbaceous understorey cover, (D) seedlings, (E) saplings, and (F) th presence of exotic species. The line shows the mean estimate, with the shading representing 50% and 90% confidence intervals. Inner and outer shading represents 50% and 90% confidence levels respectively.

We found deer density to have several plausible relationships with several structural vegetation componentsw, including negative relationships with bare ground cover (\(\beta\) = -0.04 [90% CI: -0.088, 0.005]), native woody understorey cover (\(\beta\) = -0.046 [90% CI: -0.09, -0.004]), and a positive relationship with native non-woody/herbaceous understorey cover (\(\beta\) = 0.066 [90% CI: 0.02, 0.113]). We also found strong evidence suggesting the presence of weeds is more likely at sites with higher densities of deer (\(\beta\) = 0.157 [90% CI: 0.041, 0.281], \(OR\) = 1.17 [90% CI: 1.042, 1.325]). We did not find tangible relationships between deer abundance and seedling or sapling counts.

Analysis Summary

Our studies detected deer at 148/317 cameras across Victoria and some form of deer sign (camera or transect) at 186/317 sites. Our average estimates of deer density at the 317 sites ranged from 0 to 28.3 deer per km2. Across approximately 74,570 km2 of public land in Victoria we estimate total deer abundance of the four species investigated in this study (Sambar, Fallow, Red and Hog) at 191,153 (90 % CI: 146,732 - 255,490. We also find plausible links between deer abundance and increased presence of exotic weed species.

Comparisons to other studies

Below we show a table of studies that have previously estimated deer density in Australia.

Show code
abundance_studies <- readr::read_csv("data/Australian_Deer_Studies_Density_Estimates.csv") %>%
  mutate(`Density LB` = as.numeric(stringr::str_trim(`Density Min`)), 
            `Density UB` = as.numeric(stringr::str_trim(`Density Max`)))

abundance_table_present <- abundance_studies %>%
  transmute(State, Species, Title, Author, `Study Sites`, `Year Start`, 
            `Year End`, `Density Mean`, `Density LB`, 
            `Density UB`) %>%
  arrange(State, Species, Title, `Study Sites`)

write.csv(abundance_table_present, "outputs/tables/other_studies.csv", row.names = FALSE)

abundance_table_present %>% 
  kbl("html", digits = 3) %>%
  kable_styling("striped") %>%
  collapse_rows(1:3)
State Species Title Author Study Sites Year Start Year End Density Mean Density LB Density UB
NSW Fallow Estimating deer density and abundance using spatial mark–resight models with camera trap data Bengsen, A. J., Forsyth, D. M., Ramsey, D. S. L., Amos, M., Brennan, M., Pople, A. R., Comte, S., & Crittle, T Montane Woodland/Fores in Kosciuszko NP 2017 2019 0.29 0.12 0.53
Sambar Bengsen, A. J., Forsyth, D. M., Ramsey, D. S. L., Amos, M., Brennan, M., Pople, A. R., Comte, S., & Crittle, T Montane Forest in Kosciuszko NP 2017 2019 0.73 0.33 1.35
Bengsen, A. J., Forsyth, D. M., Ramsey, D. S. L., Amos, M., Brennan, M., Pople, A. R., Comte, S., & Crittle, T Montane Woodland/Fores in Kosciuszko NP 2017 2019 2.49 1.67 3.50
Bengsen, A. J., Forsyth, D. M., Ramsey, D. S. L., Amos, M., Brennan, M., Pople, A. R., Comte, S., & Crittle, T Wetland Reserve 2017 2019 0.48 0.24 0.84
QLD Red I just want to count them! Considerations when choosing a deer population monitoring method Amos, M., G. Baxter, N. Finch, A. Lisle, and P. Murray Cressbrook Dam catchment 2010 2012 45.75 39.80 51.70
Amos, M., G. Baxter, N. Finch, A. Lisle, and P. Murray Cressbrook Dam catchment 2010 2012 26.50 23.70 29.30
Amos, M., G. Baxter, N. Finch, A. Lisle, and P. Murray Cressbrook Dam catchment 2011 2011 26.26 12.92 53.30
Rusa Estimating deer density and abundance using spatial mark–resight models with camera trap data Bengsen, A. J., Forsyth, D. M., Ramsey, D. S. L., Amos, M., Brennan, M., Pople, A. R., Comte, S., & Crittle, T Central Coastal QLD Woodland 2017 2019 10.34 7.84 13.32
Bengsen, A. J., Forsyth, D. M., Ramsey, D. S. L., Amos, M., Brennan, M., Pople, A. R., Comte, S., & Crittle, T Central Coastal QLD Woodland 2017 2019 0.68 0.21 1.77
Bengsen, A. J., Forsyth, D. M., Ramsey, D. S. L., Amos, M., Brennan, M., Pople, A. R., Comte, S., & Crittle, T Samsonvale Resivoir 2017 2019 3.11 1.76 5.07
TAS Fallow Report of state-wide census of wild fallow deer in Tasmania project: Part A: Baseline aerial survey of fallow deer population, central and north-eastern Tasmania. Lethbridge, M. R., Stead, M. G., & Wells, C. Northern Midlands 2019 2019 2.70 1.67 3.71
VIC Estimating deer density and abundance using spatial mark–resight models with camera trap data Bengsen, A. J., Forsyth, D. M., Ramsey, D. S. L., Amos, M., Brennan, M., Pople, A. R., Comte, S., & Crittle, T Cardinia Resivoir Woodland 2017 2019 2.09 1.46 2.92
Hog The influence of evolutionary history and body size on partitioning of habitat resources by mammalian herbivores in south-eastern Australia Davis, N. E., Gordon, I. R., & Coulson, G Yanicke 2003 2004 3.01 1.73 4.29
Red Estimating deer density and abundance using spatial mark–resight models with camera trap data Bengsen, A. J., Forsyth, D. M., Ramsey, D. S. L., Amos, M., Brennan, M., Pople, A. R., Comte, S., & Crittle, T Sugarloaf Resivoir Woodland 2017 2019 24.57 19.79 30.64
Bengsen, A. J., Forsyth, D. M., Ramsey, D. S. L., Amos, M., Brennan, M., Pople, A. R., Comte, S., & Crittle, T Yan Yean Woodland 2017 2019 19.76 17.58 22.17
Sambar BBRR Feral Pig Survey Project Delicknora 2023 Durand M., and Tame, A. Delicknora NA NA 1.70 1.70 1.70
BBRR Feral Pig Survey Project Wulgulmerang 2022 Durand M., and Tame, A. Wulgulmerang 2022 2022 3.29 3.29 3.29
BBRR Feral Pig Survey Project Wulgulmerang 2023 Durand M., and Tame, A. Wulgulmerang NA NA 6.40 6.40 6.40
Estimating deer density and abundance using spatial mark–resight models with camera trap data Bengsen, A. J., Forsyth, D. M., Ramsey, D. S. L., Amos, M., Brennan, M., Pople, A. R., Comte, S., & Crittle, T Cardinia Resivoir Woodland 2017 2019 11.94 8.44 16.48
Bengsen, A. J., Forsyth, D. M., Ramsey, D. S. L., Amos, M., Brennan, M., Pople, A. R., Comte, S., & Crittle, T Yan Yean Woodland 2017 2019 3.93 2.53 6.29
The application of catch–effort models to estimate the efficacy of aerial shooting operations on sambar deer Ramsey David S. L., McMaster Damien, Thomas Elaine Alpine NP - Bogong 2020 2022 0.89 0.70 1.10
Ramsey David S. L., McMaster Damien, Thomas Elaine Alpine NP - Eastern Alps 2020 2022 1.40 1.08 1.76
Ramsey David S. L., McMaster Damien, Thomas Elaine Burrowa-Pine 2020 2022 0.40 0.27 0.56
Ramsey David S. L., McMaster Damien, Thomas Elaine Coopracambra 2020 2022 1.12 0.04 0.20
Ramsey David S. L., McMaster Damien, Thomas Elaine Croajingolong NP 2020 2022 0.61 0.44 0.81
Ramsey David S. L., McMaster Damien, Thomas Elaine Errinundra NP 2020 2022 0.80 0.58 1.05
Ramsey David S. L., McMaster Damien, Thomas Elaine Mount Mitta Mitta RP 2020 2022 0.30 0.14 0.49
Ramsey David S. L., McMaster Damien, Thomas Elaine Mt Buffalo NP 2020 2022 0.78 0.64 0.94
Ramsey David S. L., McMaster Damien, Thomas Elaine Snowy River NP 2020 2022 1.08 0.86 1.28
Ramsey David S. L., McMaster Damien, Thomas Elaine Wabba WP 2020 2022 1.30 1.02 1.60
The role of wild deer in the transmission of diseases of livestock Pacioni, Ramsey et al Gembrook 2018 2020 5.90 5.70 6.10
Pacioni, Ramsey et al Kinglake 2018 2020 6.80 3.70 10.90
Pacioni, Ramsey et al Willow Grove 2018 2020 6.20 5.90 6.50

Presence-absence records

The camera traps used in this study also detected numerous other species. We can obtain a table of the presence-absences of these species using the following code:

Show code
# Get view 
pa <- tbl(con, dbplyr::in_schema("camtrap", "processed_site_substation_presence_absence")) %>%
  dplyr::filter(ProjectShortName %in% !!project_short_name) %>% # Only retrieve cam records for hog deer 2023
  dplyr::collect() %>%
  dplyr::arrange(SiteID) 

# Get pa table
pa_table <- pa %>%
  group_by(Species = common_name) %>%
  summarise(`Sites Detected` = sum(Presence)) %>%
  arrange(desc(`Sites Detected`))

write.csv(pa_table, "outputs/tables/other_species.csv", row.names = F)
Bürkner, Paul-Christian. 2017. brms: An R Package for Bayesian Multilevel Models Using Stan.” Journal of Statistical Software 80 (1): 1–28. https://doi.org/10.18637/jss.v080.i01.
———. 2018. “Advanced Bayesian Multilevel Modeling with the R Package brms.” The R Journal 10 (1): 395–411. https://doi.org/10.32614/RJ-2018-017.
———. 2021. “Bayesian Item Response Modeling in R with brms and Stan.” Journal of Statistical Software 100 (5): 1–54. https://doi.org/10.18637/jss.v100.i05.
Vehtari, Aki, Jonah Gabry, Mans Magnusson, Yuling Yao, Paul-Christian Bürkner, Topi Paananen, and Andrew Gelman. 2020. “Loo: Efficient Leave-One-Out Cross-Validation and WAIC for Bayesian Models.” https://mc-stan.org/loo/.
Vehtari, Aki, Andrew Gelman, and Jonah Gabry. 2017. “Practical Bayesian Model Evaluation Using Leave-One-Out Cross-Validation and WAIC.” Statistics and Computing 27 (5): 1413–32.

References